Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

Low-frequency wave-energy amplification in graded two-dimensional resonator arrays



    Energy amplification in square-lattice arrays of C-shaped low-frequency resonators, where the resonator radii are graded with distance, is investigated in the two-dimensional linear acoustics setting for both infinite (in one dimension) and finite arrays. Large amplifications of the incident energy are shown in certain array locations. The phenomenon is analysed using: (i) band diagrams for doubly-periodic arrays; (ii) numerical simulations for infinite and finite arrays; and (iii) eigenvalue analysis of transfer matrices operating over individual columns of the array. It is shown that the locations of the large amplifications are predicted by propagation cut-offs in the modes associated with the transfer-matrix eigenvalues. For the infinite array, the eigenvalues form a countable set, and for the low frequencies considered, only a single propagating mode exists for a given incident wave, which cuts off within the array, leading to predictive capabilities for the amplification location. For the finite array, it is shown that (in addition to a continuous spectrum of modes) multiple discrete propagating modes can be excited, with the grading generating new modes, as well as cutting others off, leading to complicated amplification patterns. The numerical simulations reveal that the largest amplifications are achieved for a single row array, with amplifications an order of magnitude smaller for the corresponding infinite array.

    This article is part of the theme issue ‘Modelling of dynamic phenomena and localization in structured media (part 1)’.

    1. Introduction

    Recently, Bennetts et al. [1] have shown that C-shaped bottom-mounted cylinders act as local low-frequency (sub-wavelength) resonators in a line array subject to plane incident waves. A spatial grading of the cylinder radii led to a structure amplifying the incident wave energy by factors up to more than 500—over 11 times greater than the amplification in an isolated cylinder—and over a broad range of wavelengths, with the position of the largest amplification in the array depending on the incident wavelength. The effect is known as rainbow trapping after Tsakmakidis et al. [2], who illustrated, in optics, how grading the material properties of a tapered waveguide could lead to spatial separation of its spectrum and trapping of the wave; subsequently this idea has been used to advantage in many areas of wave physics, e.g. acoustics [3], elastic waves [4] and plasmonics [5], among others. To illustrate this effect, following ref. [1], figure 1 shows the energy, E, normalized to unit incident energy for a graded row of 10 C-shaped resonators and an incident wavenumber close to the resonant wavenumber of the fifth resonator (counted from the left). The maximum energy amplification occurs in the fourth cylinder; it is more than 270 times greater than the incident wave energy, and 6.9 times the amplification in the isolated cylinder.

    Figure 1.

    Figure 1. Logarithm of normalized energy (unit incident energy) along a single graded row of 10 C-shaped resonators for am/a1 = (m + 8)/9, a1/W = 0.217, φ = 0.1 π, kW = 0.476 π and ψ = 0.

    While the results given in Bennetts et al. [1] were in the three-dimensional water-wave context, motivated by supplying design ideas for wave-energy-conversion structures, the results can be directly interpreted in a two-dimensional acoustics setting, which will be considered in what follows. Here, the field satisfies the Helmholtz equation (with wavenumber k) in the x = (x, y) plane and the C-shapes are sound-hard structures (i.e. homogeneous Neumann conditions are assumed on their surfaces). More precisely, the C-shapes are defined by

    for m = 1, …, M, where M = 10, (xm, 0) is the centre location of the mth C-shape, W = xm+1 − xm is the constant spacing (centre to centre), am is the radius of the mth C-shape and φ is the half-angle of its opening (identical for all C-shapes), with the openings at the left-hand end and symmetric about the line y = 0, as shown in figure 1. Correspondingly, the Helmholtz equation is imposed in R2mCm¯. Motions are forced by a plane incident wave, with velocity potential
    where ψ is the incident wave direction with respect to the positive x-axis (ψ = 0 corresponds to head-on incidence, as in figure 1). The system is completed by the Sommerfeld radiation condition as well as a condition ensuring the correct singularity at the tips of the C-shapes.

    It was shown in Bennetts et al. [1] that the rainbow-trapping effect is explained by the behaviour of local Rayleigh–Bloch waves excited in the array. The C-shape-radius grading causes the rightward-propagating Rayleigh–Bloch wave excited at the leading end of the array to slow down progressively along the array, until it reaches a turning point, at which the group velocity is zero. In the language of waveguide modes, at this point, the wave is cut off and ceases to propagate. The energy carried by the local Rayleigh–Bloch wave accumulates at the turning point, thereby generating large amplifications.

    Graded two-dimensional arrays of structures have been shown to support rainbow trapping. In an acoustic setting, Romero-Garcia et al. [6] showed that a graded array of ordinary cylinders (non-resonant closed circles; in the context introduced above, a zero opening angle, φ = 0), with lattice spacing decreasing in the direction of the incident wave, can be used to amplify acoustic pressure in certain regions. They showed this experimentally but also argued that the effect is due to the incident wave reaching a region inside the grating where its effective group velocity vanishes so that the incident energy accumulates. Band-structure calculations were used to support this argument.

    The present work is an extension of Bennetts et al. [1] to two-dimensional graded arrays on square lattices, in which the locations of the C-shapes are

    for m = 1, …, M and n = 1, …, N, where W = xm+1 − xm = yn+1 − yn. It can also be viewed as an extension of Romero-Garcia et al. [6] to resonating members with grading in terms of the resonance properties. Section 2 presents band structures of doubly-periodic arrays of C-shaped resonators, from which the effects of grading can be inferred. Section 3 considers effects of grading for stacks of diffraction gratings, where each grating contains an infinite number of identical C-shaped resonators arranged periodically in a single column. A transfer-matrix method is outlined and the grading effects are related to changes in the spectrum of the underlying (discrete) transfer operator. Section 4 considers the associated finite-array problem, in which each diffraction grating is truncated to contain only a finite number of identical C-shaped resonators. An analogous transfer-matrix method is employed, which is now based on discretization of a transfer operator associated with a continuum of directions. Section 5 ends the paper with a short summary and brief discussion of the findings.

    2. Band diagrams

    In order to understand the behaviour of waves travelling through a graded two-dimensional array of structures on a square lattice, it is enlightening to consider the periodic case—a so-called doubly-periodic array, i.e. no grading—the idea being that, as long as the grading is weak enough, the structure will approximately experience a wave field as if it were a member in a periodic array of (global) periodicity equal to the actual local periodicity.

    In the periodic case, the problem reduces to finding solutions satisfying the Bloch condition

    for all lattice vectors R = (m1W, m2W), where m1 and m2 are integers and q = (q1, q2) is the Bloch wavevector, which is part of the unknowns. This problem can be written as an eigenvalue problem with quasi-periodic boundary conditions on the unit cell ( − W/2, W/2)2. More precisely, the wavenumber in the Helmholtz equation is now the (square-root of the) eigenvalue of the Laplacian subject to (Floquet–)Bloch conditions at the sides of the unit cell, which are given by
    as well as the (unchanged) homogeneous Neumann condition on the C-shaped surface.

    As we are interested in propagating waves, only real-valued entries of q are considered. Moreover, by periodicity, it suffices to consider values of q in the so-called first Brillouin zone {qW | qW∈[ − π, π)2}. In fact, owing to the symmetry of the C-shapes about the x-axis, q2 is symmetric about the q1-axis and, by a time-reversal argument, the eigenvalue also does not change if q is swapped for −q. The full picture is thus given by {qW | qW∈[0, π)2}, which makes up the so-called irreducible Brillouin zone in this case.

    We solve the Bloch problem using the finite-element method. In order to simplify the numerical computation, we add a small thickness to the C-shaped walls. Thus, the radius in the following results refers to the inner radius (decisive for the internal resonance). Results for the fifth C-shape (radius = a5) are displayed in figure 2. Figure 2a shows the first and the second band (solid blue lines) for q2 = 0. In solid-state physics, the left-hand end, q1W = 0, is called the Γ-point, while the right-hand end, q1W = π, is referred to as the X-point. For comparison, the first band of the equivalent closed circle (φ = 0, same outer radius) is also shown (solid grey curve) as well as the so-called light line (dashed black line), which is the band of free space. As can be seen, the presence of the scatterers lowers the wavenumber in comparison to free space and the C-shaped resonators do this much more so than the closed circles. Moreover, a band gap for waves propagating in the x direction is seen for approximately kW/π∈[0.4, 0.6].

    Figure 2.

    Figure 2. (a) Dispersion curves for q2W = 0, a/W = 0.313 and φ = 0.1 π (i.e. C-shaped resonator; solid blue curves), φ = 0 (i.e. closed circles; solid grey curve), with the light line overlaid (dashed black line). (b) Dispersion surfaces for the C-shaped resonator (φ = 0.1 π) in the first Brillouin zone.

    To investigate the band gap for waves propagating in other directions, the first two dispersion surfaces for the C-shape are plotted for the full first Brillouin zone in figure 2b. (As noted above, it would have been sufficient to have the plot for the irreducible Brillouin zone; in particular, the symmetry alluded to above is confirmed by the figure.) The first band goes beyond its value at the X-point and the maximum occurs at the point (q1W, q2W) = (π, π), the so-called M-point. Therefore, the all-angle band gap is approximately kW/π∈[0.45, 0.6]. Note that the resonant wavenumber of the equivalent closed circle is kW/π ≈ 0.476, which induces the cut-off for the first band.

    As the C-shapes become smaller in radius, the related resonant wavenumber becomes larger and the cut-off moves up in the band diagram. Accordingly, the first bands of smaller C-shaped-resonator radii go to higher wavenumbers (not shown). Therefore, a gently graded array of C-shaped resonators, where the radius of the C-shapes increases away from the front of the array, is expected to allow waves to penetrate into the array up to the location at which the associated band diagram has a gap for the wavenumber of the incident wave. Note also that the slope of the dispersion curve gives the group velocity of the wave, which is the velocity at which energy propagates. As the slope of the first band is zero at the X-point, the group velocity is zero at this point, causing energy to accumulate.

    3. Stacks of infinite diffraction gratings

    The key to understanding the effects of grading on a stack of infinite diffraction gratings lies in analysing how the eigenvalues of the transfer matrices for the individual gratings evolve with the grading. In each of the M gratings, the C-shaped resonators are repeated periodically in the ±y direction, so that they form the columns of the array. Reflection and transmission of a plane incident wave by a single grating is calculated using the method given by, for example, Peter et al. [7], with the diffraction transfer matrix for a single C-shaped resonator calculated using the method of Montiel et al. [8], in which C-shaped resonators are treated as being infinitely thin (i.e. one-dimensional objects). Wave interactions between the gratings in a stack are then calculated using the method developed elsewhere [911].

    We briefly lay out the ideas in the notation of this paper. The scattered solution for a single grating can be expressed as a countable sum of propagating/decaying plane waves. As in the finite-array case introduced in §1, the distance between centres of adjacent C-shaped scatterers is denoted by W, and the incident wave, ϕinc, propagates towards the array at an angle ψ with respect to the positive x-axis, as in (1.2). The directions of the scattered waves are determined by the discrete set of real solutions ψj of

    The elements of the set {ψj:jZ} are referred to as the scattering angles. They depend only on the y-periodicity dictated by the incident wave and the spacing W, and are independent of the scatterer properties (i.e. C-shaped resonator properties).

    As indicated in (1.3), the columns (i.e. gratings) are labelled m = 1, …, M, and ordered left to right. Column m occupies the x interval xm < x < xm+, where xm ±  = xm ± W/2, and the centres of the C-shaped resonators in the column lie on the contour x = xm, noting that xm+ = xm+1− for m = 1, …, M − 1. The wave field between columns m and m + 1 is

    where a(m ± )j and b(m ± )j, respectively, are the incoming and outgoing wave amplitudes, from the left ( − ) and right ( + ) of column m.

    It is convenient to introduce matrix and vector notation at this point. Therefore, following truncation, expressions (3.2) are recast as

    using the following notation for diagonal matrices,
    and for vectors,

    The transfer matrix Pm is the matrix that relates amplitudes on either side of column m, i.e.

    The eigenvalues of the transfer matrix, {μC:μeig(Pm)}, make up the spectrum of the associated operator, and define the waves supported by the corresponding doubly-periodic array. In particular, eigenvalues on the unit circle, |μ| = 1, define Bloch wavenumbers on the bands via the relation μ = exp(i q1W), with q2 defined by the incident wave via q2 = k sinψ.

    Figure 3a shows results analogous to those of figure 1 but for a single row within a stack of infinite gratings. In this case, the maximum amplification occurs in the third column of C-shaped resonators, m = 3, as opposed to the fourth column in figure 1. The amplification is weaker in the infinite-row array (stack of gratings) than the single-row array shown in figure 1: it is ≈50 times the incident energy and 2.6 times the amplification of the isolated cylinder. Figure 3b shows results for an incident wave with angle ψ = π/4. The maximum amplification also occurs in the third column but is slightly weaker, at ≈30 times the incident energy and 1.6 times the amplification of the isolated cylinder. The amplification is almost the same in the fourth column, indicating that the wave has penetrated farther into the array before reaching the cut-off. The behaviour is consistent with the band diagrams shown in figure 2, in which, at the M-point, the dispersion surface extends beyond its level at the X-point. (Note that these diagrams are for the fifth C-shape but analogous results hold for the other ones.)

    Figure 3.

    Figure 3. (a) As in figure 1, but for a single row of a stack of infinite diffraction gratings. (b) Corresponding results for oblique incidence with ψ = π/4. (cg) Dominant eigenvalues of transfer matrices in complex plane, for ψ = 0 (red bullets) and ψ = π/4 (blue stars) operating over columns n = 1, …, 5. Unit circles are overlaid (broken curves).

    Figure 3cg shows the eigenvalues of transfer matrices Pm for m = 1, 2, 3, 4 and 5, respectively, corresponding to the dominant modes in the respective columns. Eigenvalues are shown for ψ = 0 (red dots; corresponding to figure 3a), and for ψ = π/4 (blue stars; corresponding to figure 3b). The eigenvalues come in complex-conjugate pairs, with the eigenvalues in the upper half complex plane defining rightward-propagating waves, and the eigenvalues in the lower half plane defining leftward-propagating waves. Figure 3c shows the eigenvalues on the unit circle for the transfer matrices that map over the first column, i.e. there is a single rightward- (or leftward-) propagating mode supported by the corresponding doubly-periodic array for both ψ = 0 and ψ = π/4. Figure 3dg shows that, as the radius increases along the array columns, the eigenvalues initially move around the unit circle towards −1, with the eigenvalues for ψ = π/4 trailing behind those for ψ = 0. As the eigenvalues approach −1 from above (rightward propagating) and below (leftward propagating), the group velocity of the waves reduces, and when the eigenvalues meet at −1 they jump onto the negative real axis, i.e. the associated modes cut off and decay with no energy propagation; one eigenvalue jumps inside the unit circle (rightward decaying), and the other outside (leftward decaying; beyond axes limits). The eigenvalues for ψ = 0 cut off in column 4, whereas the eigenvalues for ψ = π/4 cut off in column 5. Therefore, the energy carried by the wave modes accumulates around the third to fourth columns, generating the amplifications shown in figure 3a,b. As the eigenvalues for ψ = π/4 take an additional column to reach the cut-off, the amplification is spread more evenly over columns m = 3 to 4. The behaviour is qualitatively similar to the single-row case shown in figure 1, in which the eigenvalues that jump off the unit circle onto the negative real axis are associated with so-called Rayleigh–Bloch modes.

    4. Finite arrays

    We consider stacks of finite gratings, in which the gratings are truncated after a finite number of C-shapes (symmetrically about the x-axis). This provides an intermediate situation between the columns containing a single C-shaped resonator (as in figure 1) and an infinite grating (as in §3). The analysis can be performed in a similar fashion to the infinite gratings in §3 but the reduction to a finite number of resonators in the column complicates the situation. In particular, waves propagate in all directions for the finite grating, which means that the analogue of the transfer matrix from §3 is an operator acting on a continuum of directions. The number of discrete eigenvalues is broadly related to the number of members in the column. This directly relates to the transfer-matrix method recently presented in Bennetts et al. [12], which in turn is based on Montiel et al. [13,14]. We briefly summarize the essence of the method in what follows. A version of this method for a single resonator was used by Bennetts et al. [1] to analyse the results of figure 1.

    Consider column m, which occupies the region

    Outside of the resonators, the total wave field, ϕ, is decomposed into fields propagating/decaying rightwards and leftwards via the integral representations
    ϕ(x,y)=Γ±Am±(χ)φm±(x,y:χ)dχ+ΓBm±(χ)φm±(x,y:χ)dχ,for ±(xxm)0,4.2
    is a plane wave travelling in direction χ and normalized to the left-hand side (−) or right-hand side (+) of the column. The integration contours are defined by
    and Γ+ = Γ + π. On the real branches of Γ ± , φm ±  defines rightward-propagating waves for Γ and leftward-propagating waves for Γ+. On the complex branches, φm ±  defines plane waves that decay rightwards for Γ and leftwards for Γ+. Note that (4.2) for + is the analogue of (3.2) but for a continuum of propagation angles.

    In (4.2), Am ±  and Bm ±  represent incoming and outgoing amplitude functions, respectively. The outgoing amplitude functions are expressed in terms of the incoming amplitude functions using the scattering relations

    Here, Rm and Tm are, respectively, the reflection and transmission kernels for column m. They determine the outgoing amplitude response in direction χ due to an incoming amplitude forcing in direction ψ.

    The transfer operator is now precisely the operator mapping the amplitude functions from the left boundary of the column across to the right boundary, i.e.

    Note the close similarity of (4.6) to (3.4). However, (3.4) is a relation for a discrete number of scattering angles, while (4.6) is a relation for a continuum of directions and is independent of the incident field, as there is no imposed quasi-periodicity, i.e. no q2. Relation (4.6) becomes a finite-dimensional matrix relation after discretization (and truncation) of the complex contours Γ±.

    Figure 4ac shows spatial energy distributions for a graded array of 10 × 10 C-shaped resonators, in which each column is made up of 10 identical resonators of the type shown in figure 1, and for (a) normal incidence ψ = 0, (b) oblique incidence ψ = π/4 and (c) grazing incidence ψ = π/2. For normal incidence, the maximum energy amplification is approximately 104 times the incident energy, which is more than twice the maximum amplification in the corresponding infinite array. The maximum occurs in the fifth column and the second and ninth rows (counted from the top), i.e. for cylinders that resonate in isolation for the chosen incident wavenumber; the amplifications are 1.8 times greater than for the cylinders in isolation. Large amplifications also occur in the third column and the fourth and seventh rows, for which the energy amplifications are ≈87 times the incident wave energy and 4.6 times the amplifications in the isolated cylinders; and in the third and sixth rows and the fourth column in rows two and nine, for which the amplifications are ≈70 times the incident wave and 1.7 times the isolated cylinders. Hardly any energy reaches the sixth column (maximum amplification of about 2 occurring in rows two and nine), and virtually none in columns farther into the array. Boundary effects are visible at the upper and lower edges of the array. For both oblique and grazing incidence, the amplification patterns are broadly similar, with cut-offs after the fifth column. The energy amplifications are generally larger: more than 210 times the incident wave energy for ψ = π/4 and more than 185 for ψ = π/2. Note that, for a single row (N = 1), as in figure 1, the discrete spectrum contains a single complex-conjugate eigenvalue pair corresponding to the Rayleigh–Bloch wave [1], which cannot be excited at grazing incidence due to symmetry arguments.

    Figure 4.

    Figure 4. (a) As in figure 1 but for 10 rows (N = 10). (b) As in (a) but for ψ = π/4. (c) As in (a) but for ψ = π/2. (dh) Eigenvalues of transfer matrices corresponding to point spectrum of columns m = 1, …, 5.

    Figure 4dh shows eigenvalues of the transfer matrices Pm for m = 1, …, 5, respectively. In contrast to the eigenvalues for infinite gratings shown in figure 3cg, the eigenvalues for finite gratings are the same for normal, oblique and grazing incidence. Only the point eigenvalues are shown; for the first column (smallest resonators), there are four eigenvalues on the unit circle in the upper half plane (and four in the lower half plane), which are fractionally displaced along from the ends of the continuous spectrum (the numerical approximation of which is omitted for clarity in the plots). The eigenvalue farthest away from the continuous spectrum is labelled μ1 and the closest is μ4, as shown in figure 4d. In a similar fashion to the eigenvalues of the infinite grating shown in figure 3cg, as the resonators get larger from m = 1 to 5, the eigenvalues move around the unit circle towards −1 before jumping onto the negative real line; the locations of μ1 and μ4 for the second and fourth columns are indicated on figure 4e,g, respectively. Therefore, analogous conclusions about the cut-off of modes for the infinite gratings in §3 can be drawn for the finite array. However, the behaviour is more complicated for the finite array, because (i) there are more (point) eigenvalues and (ii) additional eigenvalues appear from the ends of the continuous spectrum as the resonators get larger, with five eigenvalues in the upper half plane for column m = 2, and nine in the upper half complex plane/on the negative real axis within the unit circle for column m = 5. This leads to non-uniform response in the different rows (even for normal incidence), and maximum amplifications appearing in different rows along the array.

    In order to gain further insight into the resonant behaviour, figure 5 shows energy amplifications for the normal-incidence wave field (figure 4a) in columns m = 4 and 5 projected onto the most strongly excited modes Mj, associated with the eigenvalues μj labelled in figure 4g,h. In the fourth column, M1 is the dominant mode, followed closely by M4, and the left-hand and middle panels of figure 5a show their respective contributions to the wave field. Mode M1 has cut-off in the fourth column, and M4 is the next mode to cut off after the fourth column. The right-hand panel of figure 5a shows the sum of all modes combined, 17, where the sum includes the contributions of the leftward-propagating/decaying modes in the point spectrum from which the dominance of M1 and M4 is clearly seen. Figure 5b shows analogous results for the fifth column of the array, in which the dominant modes are M7 and M9, with M7 having cut-off in the fifth column and M9 about to cut off. Again, it is clearly seen from these plots that the response in the column is dominated by these modes.

    Figure 5.

    Figure 5. Left-hand and middle panels show normalized logarithmic energy supported by the two dominant modes excited by a normally incident wave, as shown in figure 4a, for (a) column m = 4 and (b) column m = 5. Right-hand panels show combined contribution of all modes corresponding to discrete spectrum.

    5. Conclusion

    Two-dimensional square-lattice resonator arrays, in which the resonator radii are graded with distance from the front of the array, have been shown to generate large energy amplifications at certain locations in the array. Understanding of the amplification locations was derived by treating the array as a stack of diffraction gratings, and analysing the locations along the array at which the modes associated with eigenvalues of the transfer matrices acting over the individual gratings cut-off. It was shown that the maximum amplifications are larger for finite arrays than infinite arrays, but that the amplification patterns are more complicated for finite arrays. The complicated amplification pattern was related to excitation of multiple modes corresponding to point eigenvalues, with different cut-off locations along the array.

    Data accessibility

    Data used in this article are available from the first author.

    Authors' contributions

    All authors contributed towards the investigation, including reading and approving the manuscript.

    Competing interests

    The authors declare that they have no competing interests.


    R.V.C. thanks the UK EPSRC for their support through Programme Grant EP/L024926/1 and also acknowledges the support of the Leverhulme Trust.


    M.A.P. gratefully acknowledges the hospitality of the University of Adelaide, at which most of this paper was written during a research visit. The authors thank Fabien Montiel for providing the codes to calculate scattering by an individual C-shaped resonator and finite arrays.


    One contribution of 14 to a theme issue ‘Modelling of dynamic phenomena and localization in structured media (part 1)’.

    Published by the Royal Society. All rights reserved.