Beam path and focusing of flexural Lamb waves within phononic crystal-based acoustic lenses

We describe an analytical approach that allows calculating the trajectories of an elastic beam within a two-dimensional phononic crystal-based gradient index medium. The formalism takes into account the anisotropy along the lines of inclusions where the equi-frequency contours may depart from a circle. We then report on a numerical and experimental study of the focusing of flexural Lamb waves in gradient-indexed phononic crystals. The silicon/air heterostructures that we considered for this work features an index that is obtained through the modulation along one axis of either the diameter of the air inclusions or their spacing. In both cases, numerical and experimental results agree very well. The formalism that we have developed explains well the location, shape and size of the focus in either system.


Introduction
Focusing of elastic waves via flat phononic lenses is definitely one of the most striking phenomena that arise from the artificial periodicity of phononic crystals (PCs). After the first experimental demonstration of the phenomenon was achieved in 2008 [1], an intense activity has developed to optimize the geometrical and physical parameters of acoustical lenses in order to realize the ultimate spatial resolution [2][3][4][5][6][7]. These promising perspectives for acoustic imagery relate to the fact that PCs are systems in which one can control the wave propagation at the wavelength scale, making it possible to overcome the diffraction limit. Actually, it is not possible for a conventional lens to produce an image containing details that are finer than half the wavelength of the light or of the sound being focused. This drawback is called the 'Abbé diffraction limit', and applies whatever the waves are, electromagnetic and elastic as well. This is the consequence of the evanescent waves emerging from the object but that do not contribute to the resolution of the image because of the exponential decay of their amplitude in any lens made of a positive index material. In contrast, in lenses where negative refraction comes about, the evanescent waves have their amplitude that increases during the transmission through the medium. After emerging from the lens, their amplitude decays again down to their initial level that is reached in the image plane. However, because they can still contribute to the resolution of the image, the Abbé diffraction limit is overcome. Another consequence of the exponential decay of the evanescent waves emitted by the source is that the super-resolution applies only to near-field imaging. Actually, if the source is located at a distance from the lens that is greater than one wavelength, the decay turns out to be prohibitive for the evanescent component to reach the lens and the super-focusing effect cannot be achieved; this may constitute a severe limitation for potential applications.
At the starting point of the super-focusing effect is the negative refraction of ultrasounds that can be observed either in elastic metamaterials or in PCs. The former are constituted by a set of local resonators embedded into a matrix; the dimensions of these resonators must be much less than the acoustic wavelength in the host media, so that the effective refractive index can be defined through homogenization theories. However, for this condition to be fulfilled, the local resonators must be coupled to the matrix through a medium with very specific physical properties, making these artificial structures hard to elaborate. As regards the PCs, the negative refraction of elastic waves can be controlled through the band structure since this phenomenon is the direct consequence of the bands folding and therefore of the negative slope of some acoustical branches [8]: only elastic waves with frequencies above the first band gap can contribute to the superlensing effect. Their wavelengths are then smaller than the period of the PC lens and effective properties cannot be defined.
These two-dimensional (2D) systems are engineered with a gradual variation of their constitutive parameters (e.g., filling factors, geometry of the inclusions or material properties) along one direction. As a result, they feature a sound velocity gradient along that direction, making the focalization of an incident wave possible. Actually, when an acoustic beam propagates through a 2D GRIN PC, it encounters redirection at every virtual interface between layers, resulting in successive reorientations of the acoustic beam inside the structure. Thus, by gradually modulating the parameters of a GRIN PC, one may create a focusing trajectory for the acoustic waves [9][10][11][12][13][14][15][16][17]. In principle, this trajectory can be analytically calculated, at least for some forms of the gradient. However, deviation between the focal distance predicted by the theory and that derived from numerical simulations may be sometime significant, even in the homogenization frequency range [9][10][11]. There are several reasons for this. First, the actual 2D acoustic lenses feature discretized indices, which may be imperfectly represented by a continuous gradient. However, this is probably not the most relevant reason since small deviations of a few percent were observed for a wavelength only five times larger than the period [9][10][11]. The observed disagreement between the theory and the numerical simulations can be explained to a larger extent by the overall shape of the equi-frequency contours (EFCs). Actually, in PCs with large filling factors, the EFCs may depart from a circle, even at low frequency [9][10][11], so analyzing the trajectories in terms of an effective index may not be relevant. Instead, the effects of anisotropy are better described by considering both the group velocity and the k vector as the local parameters [21]. Based on the same idea, a Hamiltonian optics approach has been recently proposed to study light propagation in graded PCs in the short-wavelength regime [22]. However, up to now, there have only been a few theoretical works able to perform a quantitative analysis of the ray trajectory in the homogenization range [23]. This is the reason why we address this issue in this work, where we describe a method that allows analyzing the trajectory of an elastic beam propagating in a GRIN PC. Then we use this physical model to analyze the focusing properties of 2D acoustical lenses, where the velocity gradient is realized either by gradually modifying the lattice spacing or by varying the size of the air inclusions along one direction of the PC. We present then an experimental investigation of the focusing of the zero order anti-symmetric Lamb waves propagating through air hole/ silicon GRIN PCs with either of the two geometries and we compare their efficiency in focusing flexural Lamb waves at sub-wavelength.

Ray trajectory analysis
Several refractive index profiles designed for focusing or collimating an elastic wave that enters a GRIN PC have been proposed [9][10][11][12][13][14][15][16][17][18][19][20]. Deriving the trajectory y x ( ) of the waves without carrying out any approximation is generally not doable. There is, however, an index profile n y ( ) allowing for an exact determination of the acoustical rays: the transverse hyperbolic secant profile, which can be formally written as [24] where n 0 is the index along the x-axis, at the center of the lens and δ is the gradient parameter. A lens whose index features a hyperbolic secant profile is free of aberration, i.e., any ray normally incident on the lens converges to a single point on the axis, at the focal length f l , which depends only on δ through: The gradient index along y-axis obeying to equation (1) may result from the gradual variation of the filling factor. This can be achieved by modulating either the diameter of the cylindrical air inclusions or the distance separating two consecutive inclusions while keeping their diameter uniform. In the former case, the GRIN PC has a square lattice whereas the unit cell is rectangular in the latter case. Generally, as long as the wavelength is much larger than the size of the unit cell, the medium can be considered as homogeneous and the group velocity can be simply derived from the derivative of the corresponding dispersion curve at any location in the GRIN PC. This allows in turn defining an effective index n eff as being the average of the index along ΓX and along ΓM directions [4,[9][10][11][12]25].
We used a finite element method to compute the dispersion and the EFCs of A 0 Lamb mode in the case of periodic structures made of air cylinders into silicon plates having a thickness h. We first considered a rectangular unit cell whose dimensions were a along the x-axis and b along the y-axis. These axes were respectively parallel to the crystallographic directions <100> and <010> of silicon 3 . A cylindrical hole with radius r was in the center of the unit cell. To achieve an effective index with a hyperbolic secant profile, we have imposed length b to gradually vary along the y-axis, whereas both parameters a and r kept as constant values. This was equivalent to varying both the aspect ratio b a and the filling factor along the y-axis. The result for frequencies along the first branch is shown in figure 1 for different values of the ratio b a in between 1 and 2.09.
In all the calculations, the size of the unit cell along the x-axis and the thickness of the plate were a = 100 μm and h = 110 μm, respectively. The radius of the holes was kept constant to r = 40 μm. However, one should notice that, strictly speaking, GRIN PCs are not 2D phononic crystals since they do not exhibit exact periodicity along the y-axis. Nevertheless, the dispersion properties can be computed by considering as many reduced Brillouin zones as the number of discrete values taken by the parameter b. All these reduced Brillouin zones extend over that are different according to the value of b (see the inset in figure 1).
Then we considered the squared array, which was deduced from the preceding case by setting a = b. In that case, ΓX and XM have identical lengths in the reciprocal space and one obtains a dispersion curve for each value of f (figure 2). The hyperbolic secant profile was achieved by gradually varying the radius r of the holes, or equivalently, the filling factor f. Note that we fixed the different geometrical parameters in such a way that they allow for further comparisons between the two symmetries and that they remain compatible with elaboration and experimental constraints. Actually, in the calculation, we used the parameters of the lenses that we have elaborated for the experimental part of this work, namely n 0 = 1.32 and δ = 0.088 a −1 that turned out to represent the best compromise. However, anisotropy along some rows of inclusions may occur and consequently the EFCs may depart from being circular along these rows. In that case, equation (1) fails to describe the actual index profile since it does not account for this anisotropy. In the present case, the anisotropy remains relatively weak. Indeed, in the square lattice, the anisotropy coefficient in the direction θ between k and the x-axis, defined as the ratio η θ , is maximum for θ π = 4 where, depending on the filling ratio, it takes a value between −0.03 (r = 0.2a) and 0.07 (r = 0.4a). As for the rectangular lattice, the anisotropy coefficient has a maximum at an angle that depends on the filling ratio but it always remains below 0.09. Still, as shown below, equation (1) is not factored in the calculation of the ray trajectories, it was solely used to design both lenses.  (a) Ray trajectory in a simplified GRIN PC made of five layers. P 1 , P 2 and P 3 indicate three different positions along the ray trajectory and φ is the angle between the tangent to the ray trajectory and the x-axis; (b) the k vector gradually tilts from being horizontal at P 1 to a direction with an angle θ, which depends on x, at P 3 . The tangential component k x is conserved across the interface, as shown by the location dependent EFCs; (c) k vector in a typical EFC with group velocity v g normal to the EFC. The above representation is general and applies to any other ray trajectory in the lens.
To calculate the ray trajectory within the GRIN PC of an elastic beam normally incident on the lens, we first have defined φ as being the angle between the group velocity v g (tangent to the trajectory) and the x-axis (figures 3(a) and (c)). The filling factor varying from one horizontal layer to the next means one must compute the EFC in each layer, as depicted in figure 3(b). Because of Snell's law, which states that the component k x is conserved across the interface between two consecutive layers, the initial k vector tilts gradually as the wave propagates in the medium, from the horizontal direction to a maximum angle at the mid layer. In the general case, the EFCs are not circular [10,11] and the k vector makes an angle θ with respect to the x-axis, which is not equal to φ (figure 3(c)).
To account for this anisotropy, we state that the modulus of the k vector along a given row of inclusions located at position y, may be written as: X where ΓΧ k y ( ) is the k vector along ΓX at position y, and θ ( ) p y, is a function whose value deviates all the more from unity as the EFC departs from being circular.
The components of the wave vector k along the x-and y-axes are given by: x y To account for the effect of the local anisotropy, one must relate angles φ and θ. To this end, it is sufficient to observe that the direction vector to the tangent to the EFC is θ θ , is normal to this tangent and points in the direction given by a vector whose components are − In the direct space, v g is tangent to the trajectory and therefore: In the long wavelength limit, the value at x = 0 of the x-component of the wave vector varies continuously along the y-axis according to equation (1). Moreover, as a consequence of Snell's law, the component k x keeps a constant value all along the acoustical ray that originates at position (x = 0, y), which translates into: x x 0 Finally, we used an iterative procedure to derive the ray path from equations (3)-(7). Actually, for a given value of y and at each position along the x-axis, the angle θ was derived from equation (4) after setting k y ( ) x 0 to an initial value, allowing in turn computing the slope φ using equation (5), from which the trajectory y(x) is obtained by numerically integrating equation (6).
To complete this procedure for both symmetries that we have considered, one must make assumptions on the form of the function θ ( ) p y, appearing in equation (3).

Square lattice
In order to reflect well the anisotropy of the square lattice, the function θ p y ( , ) must account both for any departure from the circular shape of the EFC and for the four-fold axis characterizing this class of symmetry. Introducing the dimensionless coefficient α = in the definition of θ p y ( , ) allows satisfying the first requirement. On the other hand, we have introduced a cosine function in the definition of θ p y ( , ) in order to fulfill the condition of periodicity of the wave vector, as θ scans all angles within π − [0 2 ] in the first Brillouin zone. The analytical form of θ k y ( , ) that closest matches the EFCs, whatever the filling factor is within the interval − [0.07 0.5], corresponding to the actual samples that we have investigated (see below), is: X The EFCs at 5 MHz for several values of the filling factor computed using a finite element method (full lines) or derived from equation (8) (indicated by circular markers) are displayed in figure 4. The agreement between both is excellent. We have also computed the EFCs using both methods and setting the frequency to different values in between 3 and 13 MHz, i.e. close to the lowest frequency of the incomplete band gap at point X in the first Brillouin zone (see r = 0.4a in figure 2): whatever the frequency and the filling ratio are, the agreement is equally good. Deriving the trajectory y x ( ) within the GRIN PC is then straightforward. Indeed, one must simply calculate both the derivatives of k x and k y with respect to θ, using equation (4) and equation (8) and then integrate their ratio according to equation (5) and equation (6).

Rectangular lattice
As with the preceding case, one can account for the anisotropy of the rectangular lattice by introducing in the definition of θ X As can be seen from figure 5, whatever the value of the aspect ratio b a is in the sample, or equivalently the filling factor in the range [0.24-0.5] that we have investigated, the fits to the EFCs were as good as for the square lattice. Indeed, at 5 MHz, the largest deviation was less than 2%; it appears along ΓM for the aspect ratio

Numerical results
We used a finite element method (FEM) to calculate the out-of-plane component of the displacements field associated to the zero order flexural Lamb mode propagating in the GRIN PCs. a being the period along the x axis, the systems were designed to have a width along the y axis of 17a for the square lattice and 16.5a for the rectangular lattice. Both samples had a length of 120a. In addition, to best mimic the experimental situation and to avoid unwanted reflections on the boundaries of the phononic plate, the structured part of the sample was surrounded by a large area of homogeneous silicon, free of air inclusions. A line source vibrating at 5 MHz was applied in front of the GRIN PCs. We compare in figure 6 the maximum amplitude of the outof-plane displacement u z in the GRIN PC with a square lattice (figure 6(a)) and with a rectangular lattice (figure 6(b)). In these figures, the dash lines are for the ray trajectories derived from the formalism established in the preceding section.
For the GRIN PC with the square lattice, FEM simulations predict the normal displacement to be maximum at the distance 34.5a from the origin, almost twice the distance of 18a derived from equation (2). Moreover, the amplitude along the x-axis is more than 90% of its maximum value in between x = 22a to x = 40a ( figure 6(a)). These spherical aberrations that are not predicted by a simple theory of rays in a GRIN lens with a hyperbolic secant profile, are the consequence of the non-circular shape of the EFC. Actually, in this case θ and φ are not equal, all the rays do not converge on a single point, and equation (2) is not relevant to calculate the focus distance anymore. Consequently, the elastic energy spreads out along the x-axis and the focal spot covers an area larger than expected.
As expected, the second maximum appears at x = 103a, about three times the focal length of 35a. The sound ray paths displayed in figure 6 explain the large spreading out of the elastic energy along the x-axis. Indeed, the aberrations increase as the waves propagate and the maximum amplitude on this second focus spot is only ∼84% of the value measured on the first one. At the same time, the normal displacement has a smoother profile along the x-axis. Similar features are observed for the spatial distribution of displacement in the GRIN PC with rectangular symmetry ( figure 6(b)). As in the previous case, the first focal point arises around 33.5a, far from the theoretical value derived from the effective refractive index and equation (2) (about 18a). However, this geometry leads to less spherical aberration, as can be seen from the ray paths drawn in figure 6(b), and the elastic energy spreads out over a smaller distance along the x-axis: the area where the amplitude is 0.95 times the maximum value extends over 5a in the case of the rectangular lattice, against ∼8a for the square lattice. This also explains the lateral profile along the y-axis that turns out to be sharper with the rectangular symmetry than it is with the square symmetry (see below). As expected, a second focus arises at x = 92a, about three times the focal length. The maximum amplitude of this second focus is ∼0.80 times the maximum amplitude at the first focus, corresponding to a lessening similar to what is observed with the square lattice.

Experimental results
In order to allow for relevant comparisons with the theory, two GRIN PCs were elaborated in a 110 μm thick silicon plate according to the designs already described in the previous sections. We used a non-contact laser based experimental technique [5] to excite the sample into vibration and to measure the out-of-plane component of the displacement field at any location in the GRIN PC. Ultra-short light pulses issued from a frequency-doubled (532 nm) Nd:YAG were focused onto the sample after they have passed through an amplitude mask and an imaging system. As a result, a series of fringes alternately bright and dark were produced and because of the photoelastic processes, elastic waves were excited in turn. This technique allows one to finely select any k vector in the first Brillouin zone by tuning the spacing of the light fringes-or equivalently the wavelength of the elastic waves-with the imaging system. In all the experiments described in this article, the excitation zone was located a few millimeters ahead of the PC itself, in a uniform region of the sample free from any air inclusion, and we fixed the carrier frequency at 5 MHz with a spectral width of about 1 MHz. Because of the large number of fringes in the excitation spot and their length (∼2.5 mm), both the direction and the magnitude of the excited k vector were accurately defined and the elastic waves can be considered as plane waves. We recorded the time dependence of the surface displacements at any point of the sample, inside the GRIN PC using a Michelson interferometer in which the light source was a He-Ne laser. One beam of the interferometer was focused on the sample acting as one of the mirrors of the interferometer to a spot size of~5 μm, whereas the reference beam was reflected by an actively stabilized mirror. A high-speed photodiode collected the interference pattern that was then digitized at 500 MS s −1 by a digital oscilloscope. We obtained a very good S/N ratio after averaging a few hundreds of scans. The microscope and the sample were both mounted on translation stages in such a way that the probe beam could be scanned across the sample with a precision of about 1 μm. This noncontact technique allowed us to record the out-of-plane component at any location on the surface of the sample with an accuracy of a few pm and to study in detail the focalization of the acoustic waves into the GRIN PC.
We measured the wave propagation in the middle area of our GRIN PCs. We show in figure 7 the displacement field recorded at three different times after the square symmetry sample had been excited into vibration. A first focus is observed to be centered at a focal length about x = 32.5a ( figure 7(a)), which is in good agreement both with the numerical results (x = 34.5a) and the ray trajectory analysis presented in section 2. Behind the focus spot, the elastic beam is expected first to be divergent within the waveguide, as it is shown in figure 7(b), and then to re-focus on a second point located at three times the focal distance of the lens at 97.5a (figure 7(c)). This is what is indeed observed, in good agreement with the numerical results. We measured the maximum amplitude of the second focus to be about 0.68 times of the first focal spot, which is a little less than the previous numerical case.
Very similar behaviors were observed with the experimental sample featuring the rectangular symmetry (figure 8). We measured in this latter case a focal distance of ∼28a, slightly less than the focal distance of 32.5a measured in the previous sample. As expected, a second focus is formed at x = 90a with an amplitude of about 0.63 times the maximum value. This is in good agreement with the corresponding numerical simulations.
The gain factor, defined as being the ratio of the maximum displacement at the focus to the amplitude of the elastic wave measured close to the excitation area of the lens, allows for quantitative comparisons between both structures. We found a gain factor of 3.5 (corresponding to a maximum displacement of 21 pm at the focus) with the square lattice against 3.2 (a maximum displacement of 16 pm) when using the lens with the rectangular lattice. These values agree fairly well with the ones derived from the numerical simulations: 3.9 in the former case and 3.6 in the latter case. Beside the smaller gain factor obtained with the rectangular lattice, this system features a greater uniformity in the distribution of elastic energy inside the first  focus. Actually, we measured the average of the normal amplitude in the area wherein it is more than 95% of the maximum value to be ∼0.99 times of the maximum value for the rectangular lattice against ∼0.97 for the square lattice.
We show in figures 9(a) and (b) the measured maxima of the out-of-plane displacements for the square symmetry and for the rectangular symmetry respectively. Whatever the symmetry is, the amplitude is not homogeneous in between the two focal points. This is consistent with our numerical results (see figure 6); above, we ascribed this feature to the spherical aberrations at the second focus. We must also notice that the amplitude recorded within an area about 20a long along the x-axis before the second focus is slightly less than expected from the simulations (see figure 6): we measured a mean amplitude around 0.4 (in normalized units) instead of 0.5 if we refer to the simulations. Several reasons may explain this small discrepancy, including the slight coupling of the symmetric mode [10] or tiny irregularities in shape, sizes and spacing of the holes, which may induce the non-coherent diffusion of the elastic waves onto the air inclusions.
As regards the focalization efficiency of the acoustic lenses, it can be evaluated through the profiles of the focal spot along both the x-and y-axes. In the left panel of figure 10, we show as full lines the normalized longitudinal profiles achieved with the square lattice ( figure 10(a)) and with the rectangular lattice ( figure 10(b)). We show the corresponding profiles along the y-axis in figures 10(c) and (d) respectively. From these data, the full width at half maximum (FWHM) was measured to be 7λ along the x-axis. The experimental transverse size of the spot was 0.71λ for the square lattice and 0.64λ for the rectangular lattice. These values are in good agreement with the ones derived from the numerical simulations: 0.84λ in the former case and 0.75λ in the latter case, quite close to the Abbé limit of 0.5λ.
In addition, it should be noted that the formalism developed in section 2 also allows for a quantitative analysis and efficiently predicting the waist of the beam. To show that, we have drawn the ray trajectories for more than 3000 initial positions evenly distributed along the y-direction in between ±8a and we have computed the number of rays intersecting a segment of a given length, sliding along a line parallel to the y-axis at position x = 30a. This 'linear density of ray' against the position of the segment is drawn as a blue line in figures 10(c) and (d). The agreement with the experimental results and the FEM calculations is very good, except for the peaks centered at y = ±4a that are not predicted by the model.
It is important to understand the physical reasons why the lateral profile along the y-axis is sharper with the rectangular symmetry than it is with the square symmetry. Actually, the origin of a better FWHM clearly lies in the smaller aberrations along the y-axis in the former case than the ones founded in the GRIN PC with a square lattice. This follows from the dependence against θ of the anisotropy coefficient η θ , which is very different according to the symmetry. Indeed, for both symmetries and for any ray trajectory, the wave vector k lies along a direction that makes an angle θ with respect to the x-axis ranging between 0°and ∼35°(see figures 6 and 9). In this range, for the rectangular symmetry, very little η θ ( ) depends on the aspect ratio b a: as long as this ratio takes a value of no more than ∼1.5, η θ ( ) varies almost linearly from 0 at θ = 0°to a value comprised in between 0.05 (for b a = 1.37) and 0.06 (for b a = 1) at θ = 35°. This is in contrast to the situation encountered in the lens with the square symmetry. In that case, η θ ( ) varies quasi-linearly against θ as well, but with a mean slope, which is either positive or negative and strongly depends on the aspect ratio. As long as the ratio r a keeps values around ∼0.3 or less, the anisotropy of the medium is mainly that of crystalline silicon, whereas the anisotropy of the effective medium dominates for larger values of the aspect ratio. Consequently, a greater spreading out of the rays, and a broader profile in turn, occur with the square symmetry.

Conclusions
We have shown in this work that one can accurately determine the focal length, the size of the focal spot, and the displacement distribution within a GRIN PC when accounting for the overall shape of the EFCs or equivalently, for the local anisotropy if any, within a row of inclusions. Being based on a geometrical approach, the ray analysis we have presented is phenomenological in that it gives a full account of the observation although it is not derived from a theory, however, it remains valid whatever the polarization of the waves is. Whereas the paraxial ray equation [24] is well suited to accurately determine the focusing properties of a GRIN lens only when the EFCs are circular [9], this ray analysis is more general and allows accounting for both the position of the focus on the x-axis and for the extension along the y-axis. It should also be noted that only real k vectors are considered in the formalism described in section 2 and hence evanescent waves are not involved in the focusing processes.
From the experimental side, we have demonstrated the focusing of zero-order Lamb mode at 5 MHz in GRIN PCs featuring two different designs. These heterostructures are free from curved surfaces, compact, and therefore they can be integrated easily with other phononic devices. In both systems, we found very good agreement between the numerical simulations and the experimental results. In particular, we have shown that the focusing over a spot with the lateral dimension close to the Abbé limit are easily obtained with the acoustic lens with the rectangular symmetry. The anisotropy being responsible for the spreading out at the focus, one must recognize that, on average, the ultrasound beam shall be subject along the path, to less anisotropy with the rectangular symmetry than it is with the square symmetry. However, larger vibration amplitudes are obtained with the heterostructure with the square symmetry. As predicted by the numerical simulations, a second focus point was actually observed. For both of the systems we studied, the vibration amplitudes at this second focal point were more than half the vibration amplitude at the first focus.