## Abstract

We address the problem of scattering of flexural waves obeying the biharmonic equation by a stack of a finite number of gratings. We express the solution of the scattering problem for a single grating in terms of reflection and transmission matrices, incorporating the effects of both propagating and evanescent incident waves. The plane wave expansion coefficients above and below the grating are linked to multipole coefficients within the grating using the grating sums and the Rayleigh identities. We derive the recurrence procedure giving the reflection and transmission matrices of the stack in terms of those of individual layers. Trapped waves between a pair of gratings are investigated.

## 1. Introduction

The topic we address here is the scattering of flexural waves by sets of periodic scatterers or gratings. This work is complementary to our recent investigations of the band structure for doubly periodic problems (Movchan *et al.* 2007; McPhedran *et al.* 2009). There are a number of reasons why it is important to study scattering problems involving a finite number of gratings. First, in any practical situation, the assemblage of scatterers will be finite, so it is essential to know how effects such as filtering scale with the size of the system. Second, the study of wave problems for finite systems of gratings permits access to the link between scattering matrices and properties of bases of Bloch functions (see, for example, Platts *et al.*2002, 2003; Botten *et al.* 2005). This link may be exploited to deduce scattering matrices for the half space filled with grating layers or for randomized arrangements of stacks in which the gratings vary from layer to layer. A third reason is that finite sets of layers may be used to demonstrate properties such as negative refraction or autocollimation, as has been done for photonic crystals (Botten *et al.* 2005), and may also be exploited to advantage for any other type of wave that has unusual diffractive effects.

The mathematical model we use here is based on the results of analysis of a canonical spectral problem for the biharmonic operator. Specifically, we develop a rigorous scattering model describing the interaction of plane waves with a grating composed of circular inclusions or fixed point scatterers. The solution relies on the expression of the flexural wave function in terms of both sets of plane waves and multipole expansions composed of cylindrical harmonics satisfying either Helmholtz or modified Helmholtz equations. The multipole expansions about individual cylinders are made consistent with quasi-periodicity using the Rayleigh identities (Platts *et al.* 2002; Botten *et al.* 2005), and reconstruction equations are derived to transform from the multipole coefficients back to plane wave coefficients. In this way, reflection and transmission matrices can be established for any layer in a stack. We note that the reflection and transmission matrices used here take into account evanescent input waves, as well as the effects of all propagating incident waves.

Using the reflection and transmission matrices for a single grating, we then apply a recurrence procedure to obtain the reflection and transmission matrices for a finite stack. The recurrence relations take into account coupling between gratings by both propagating and evanescent waves. The structure of the recurrence relations is well adapted to the analysis of flexural wave modes trapped between gratings. These trapped modes have been investigated previously by Evans & Porter (2007). The analysis presented here includes the detailed derivation of the conservation of energy and reciprocity relations, which may be used in the evaluation of the accuracy of the numerical computations.

The structure of the paper is as follows. The governing equations, together with the plane wave and multipole representations of solutions, are discussed in §§2–4. Quasi-periodic sums over multipoles and the Rayleigh identities employing them are discussed in §5. Furthermore, in §6, we derive the reconstruction equations, which link the multipole and plane wave coefficients. Section 7 presents the general form of the conservation of energy and reciprocity relationships. The particular case of pinned inclusions of zero radius and the properties of a single grating consisting of such inclusions are discussed in §§8 and 9. Recurrence relations for stacks of finite number of gratings are derived in §10. Section 11 contains a discussion of properties of a pair of gratings consisting of pinned inclusions of zero radius and, in particular, of trapped modes that can exist between the gratings.

## 2. Governing equations

Let *W* be a solution of the scattering problem for the biharmonic operator
2.1
where , *h* denotes the plate thickness, *ρ* is the mass density, *D*=*E**h*^{3}/(12(1−ν^{2})) is the flexural rigidity of the plate, *E* is the Young’s modulus, ν is the Poisson ratio and *ω* is the angular frequency. The scattering system will be composed of circular cylindrical objects of radius *a*, although the numerical results will be confined to the limit case of *a*→0, with the Dirichlet (clamping) conditions being set on the circular contours:
2.2
The solution of (2.1) admits the representation
2.3
where
2.4
and hence *W* is the superposition of ‘waves’ with real and imaginary wavenumbers.

The grating is aligned along the *x*-axis, and the structure will be built by adding more parallel gratings and analysed using the recurrence procedure. The frame of reference is shown in figure 1.

The incident field is the plane wave of any of the following two types.

Propagating (evanescent) solution of the Helmholtz equation: 2.5 where

*χ*_{0}is real and positive for a propagating solution, and it is pure imaginary (with the positive imaginary part corresponding to a reduction of |*W*_{i,H}| with a decrease in*y*) for an evanescent solution. In terms of an angle of incidence*θ*_{i}, .Evanescent solution of the modified Helmholtz equation: 2.6 where

The total field is represented as
2.7
where *W*_{i} is the incident field of one of the two types mentioned earlier, and *W*_{s,H} and *W*_{s,M} are the scattered fields solving either the Helmholtz or the modified Helmholtz equations. The solution *W* also satisfies the Bloch quasi-periodicity conditions along the *x*-axis, as defined by the incident field.

## 3. Plane wave expansions

Let *γ*_{+} and *γ*_{−} be straight lines parallel to the *x*-axis in the upper and lower half-planes, respectively, placed at an arbitrary distance from the grating. On these lines, we shall use the plane wave expansions of *W*_{s,H} and *W*_{s,M} as follows.

On

*γ*_{+}, 3.1 where the infinite range of*p*can be divided into a finite set of propagating waves with*χ*_{p}real and an infinite set of evanescent waves with*χ*_{p}imaginary with positive imaginary part. Similarly, for the evanescent waves corresponding to the modified Helmholtz equation, 3.2On

*γ*_{−}, 3.3 and 3.4

## 4. Bessel function expansions

In the vicinity of the circular boundary *C*_{0}, we use the Bessel function expansions
4.1
and
4.2
where *A*_{n}, *B*_{n}, *E*_{n} and *F*_{n} are unknown multipole coefficients. Owing to the boundary conditions (2.2), these coefficients are related as follows:
4.3
and
4.4

The following Bessel function expansions from Abramowitz & Stegun (1965) will be needed for the incident fields of the two types shown in §2 (see equations (2.5) and (2.6)):

4.5 where

*χ*_{0}is real and positive for a propagating wave, and*χ*_{0}is imaginary with positive imaginary part for an evanescent wave.4.6 where is imaginary with positive imaginary part.

## 5. Grating sums and Rayleigh identities

The above expansions around the central cylinder *x*=0 are related to similar expansions near other cylinders *x*=*d**n* in the periodic grating by the quasi-periodic Bloch factor , where *d* is the grating period. When forming quasi-periodic Green’s functions for the grating problem, we thus encounter sums of the form:
5.1
5.2
and
5.3
5.4
In the case of *l*=0, we shall need the following result:
5.5
where the right-hand sum is a cubically convergent sum. Expressions for grating sums over Hankel functions were given by Twersky (1961), and the sums of the *K*-type are exponentially convergent and can be evaluated by direct summation.

The Rayleigh identities are derived as in Movchan *et al.* (2007) with the use of the appropriate grating lattice sums in place of the doubly periodic lattice sums (Chin *et al.* 1994). We have, with denoting the amplitude of the incident wave (2.6) and denoting the amplitude of the incident wave (2.5),
5.6
and
5.7

These identities are uncoupled for solutions of the Helmholtz and the modified Helmholtz equations. The necessary coupling between the two types of waves is provided by the boundary conditions (4.3) and (4.4). The result is an inhomogeneous system of algebraic equations with respect to the multipole coefficients *E*_{l} and *F*_{l}:
5.8
and
5.9

## 6. Reconstruction equations

The system of equations (5.8) and (5.9) can be truncated and solved for the multipole coefficients *E*_{l} and *F*_{l}. However, many quantities of physical interest require the knowledge of the plane wave coefficients above and below the grating or the stack of gratings. These can be readily evaluated once the multipole coefficients are known using the reconstruction equations, which we will now derive.

The method for constructing the reconstruction equations is based on the use of Green’s theorem applied to a test function and the part of the total field corresponding to either the Helmholtz or the modified Helmholtz equation.

For the *Helmholtz part* of the field, the test function is
6.1
where
We note that *V*_{m,H} has the opposite quasi-periodicity property to the field *W*_{H} from equation (2.3). Thus, if we apply Green’s theorem to a rectangular cell, with width *d* equal to the period of the grating and an arbitrary finite extent along the *y*-axis, the vertical sides of the rectangle together do not contribute to the line integral
Let us next consider the integral along *γ*_{+}, the segment of length *d* on the line *y*=const.>*a* (figure 1*a*):
6.2
The result is independent of *y* in the two cases μ_{m}=−*χ*_{m} and μ_{m}=*χ*_{m}.

Similar remarks apply to the integral over the lower contour *γ*_{−}, the segment on the line *y*=const.<−*a*.

The remaining contribution to the integral comes from the boundary *C*_{0} of the circular cylinder of radius *a*. While evaluating this integral, we use the expansion (4.1) together with the following expansion of the test function:
We shall also use the Wronskian formula (Abramowitz & Stegun 1965)
6.3

With the calculations done, the final results for the reflection and transmission coefficients, *R*_{m} and *T*_{m}, are
6.4
and
6.5

For the *modified Helmholtz part* of the field, the test function is
6.6
where
Its expansion in terms of the modified Bessel functions is
We proceed as before, with two possible options of , and the required Wronskian formula (Abramowitz & Stegun 1965) is
6.7
The reconstruction equations for the reflection and transmission coefficients, and , of the modified Helmholtz equation take the form
6.8
and
6.9

## 7. Conservation of energy and reciprocity relations

We now consider conservation of energy and reciprocity relations, which govern the scattering by a single grating. We stress that, while our interest here is in the case of gratings having cylindrical scatterers, the relations we derive will hold irrespective of the shape of the scatterer. We derive these relationships for both propagating and evanescent incident waves and also for both types of reflected and transmitted waves.

We use the following integral identity pertaining to two solutions *u* and *v* of a problem of time-harmonic flexural vibrations of a thin plate:
7.1
We apply this relationship to a unit cell of the grating, with ∂*Ω* comprising two horizontal segments (*γ*_{+} and *γ*_{−}), two vertical segments (*γ*_{l} and *γ*_{r}), as well as the boundary *C*_{0} of a scatterer. In the calculations that follow, *u* and *v* will be chosen to have opposite quasi-periodicity factors, so that the contributions from *γ*_{l} and *γ*_{r} will cancel. We also assume that *u* and *v* satisfy the boundary conditions (2.2) and thus the boundary integral along *C*_{0} is equal to zero. For such fields, the above integral relation reduces to
7.2
where Helmholtz and modified Helmholtz terms appear separately.

We first apply equation (7.2) to *W*=*W*_{H}+*W*_{M} and its complex conjugate and use the plane wave expansions (2.5), (2.6) and (3.1)–(3.4) to derive the following conservation of energy relation:
7.3
Here *δ*_{p} denotes the amplitude of the Helmholtz type incident wave, whereas is the amplitude of the modified Helmholtz type incident wave, both in the direction specified by *α*_{p}=*α*_{0}+2*π**p*/*d*. The set *Ω*_{H} contains Helmholtz waves of propagating type, whereas denotes Helmholtz waves of evanescent type. The set comprises all integers *p*, since all modified Helmholtz waves are evanescent. This energy relationship agrees with that given by Evans & Porter (2007) (see their formula (4.18)) in the case of a propagating incident wave.

We now consider the reciprocity relations applying to the case of the return of the reflected order *p*_{*}, as shown in figure 1*b*. These relations are derived using equation (7.2), with *u* being the original field corresponding to *α*_{0}, and *v* is the returned field corresponding to *α*_{0}′=−*α*_{0}−2*π**p*_{*}/*d*. We use the plane wave expansions as mentioned earlier, which now incorporate the following expressions for direction sines and cosines
7.4
We thus find that the integral along *γ*_{−} is zero and the integral along *γ*_{+} gives a set of reciprocity constraints for the reflection coefficients, depending on whether the incident wave in the original problem is a propagating or evanescent Helmholtz wave or (evanescent) modified Helmholtz wave:
7.5
7.6
7.7

We now consider the reciprocity relations for the return of the transmitted order *p*_{*}, as in figure 1*c*. Applying equation (7.2) to the original and returned fields (corresponding to *α*_{0} and respectively), then gives
7.8
7.9
and
7.10

## 8. Zero radius limit

In the zero radius limit of the algebraic system (5.8) and (5.9), it may be shown that all multipole coefficients, except for *n*=0, are equal to zero. For small *a*, this statement is valid to leading order approximation, in which the leading multipole coefficients *E*_{0} and *F*_{0} satisfy two equations
8.1
and
8.2
The coefficient near (*β**a*)^{−1} in equation (8.2) must be zero, which implies leading order
8.3
Substituting equation (8.3) into (8.1), we get
8.4
Taking into account that the term in square brackets in equation (8.4) is 1+*O*(*β*^{2}*a*^{2}), we arrive at the result
8.5

According to equations (6.4) and (6.5), the reflection and transmission amplitudes for waves of the Helmholtz type are then 8.6 The amplitudes for waves of the modified Helmholtz type are 8.7 according to equations (6.8), (6.9) and (8.3).

For the discussion of grating stacks, which follows, we will need to characterize the reflection and transmission properties of a single grating for the additional case of a wave incident from below rather than from above. If the value of *α*_{0} is the same for the wave coming from below, we can show from the Rayleigh identities (8.1) and (8.2) that the values of *E*_{0} and *F*_{0} are the same for the two problems. We then need to adapt the argument of §6 for the case when the incident wave is in the region *y*<0. The results for the reflection and transmission coefficients, using the superscript ‘^{−}’ to denote the incidence from below, are:
8.8
8.9

## 9. Properties of a single grating

In figure 2, we show the normalized energy reflected |*R*_{0}|^{2} as a function of *β* for a single grating of pinned cylinders in the zero radius limit. Results are given for both normal incidence and an angle of incidence of 30^{°}. The normal incidence curve shows a single Wood anomaly (McPhedran & Waterworth 1972), which is associated with the Rayleigh value *β*=2*π* (this is the value at which the orders *p*=±1 become propagating). Before this is the reflection peak at which |*R*_{0}|=1, employed by Evans & Porter (2007) in their construction of trapped waves for a pair of gratings. The curve for an angle of incidence of 30^{°} has two Wood anomalies (McPhedran & Waterworth 1972), associated with *β*=4*π*/3*d* (*p*=−1) and *β*=8*π*/3*d* (*p*=−2). Only the first of these is preceded by a point of total reflectance. Note that the Rayleigh values are associated with zero reflectance points on the curves.

These properties can be understood on the basis of the lattice sum formula (5.5) and the formula for *R*_{0} in the zero radius limit, assuming a propagating input wave:
9.1
As *β* increases, we reach successive singularities of the lattice sums at each Rayleigh value, which occur first in the imaginary part and then in the real part (figure 3). These singularities in the denominator of equation (9.1) ensure that *R*_{0} is zero at Rayleigh values. Zeros of transmittance *T*_{0} require that *R*_{0}=−1, which in turn leads to the requirement that as *β* approaches the lowest Rayleigh value, the following identity holds:
9.2
Now, for small *β*, the sum on the left-hand side of equation (9.2) is much larger than the sum on the right-hand side owing to the element , while the latter diverges as *β* approaches the lowest Rayleigh value, so the two must satisfy equation (9.2) in this *β* range.

Consider next the long wavelength limit, in which *β* tends to zero, with *α*_{0} also tending to zero in such a way as to keep the incident wave of propagating Helmholtz type. Then in the denominator of equation (9.1), the dominant terms are 1/*χ*_{0} and , leading to the estimate
9.3
This gives *R*_{0}≃(−1+i)/2, for normal incidence and , for an angle of incidence of 30^{°}.

It is interesting to compare the reflectance curves in figure 2 with those in figure 4, for the case of an incident wave of modified Helmholtz type. The latter case is essentially a rescaled (diminished) version of the former, with the denominator (larger for off-axis incidence) replacing . Evanescent incident waves cannot be generated in an unbounded region, but can occur in a finite region between interacting gratings, as we now proceed to show.

## 10. Recurrence relations for grating stacks

We consider a stack of *N*_{g} gratings with the common period *d* along the *x*-axis. The recurrence procedure adopted here for the evaluation of the transmission and reflection matrices of the stack of gratings corresponds formally to that used in Platts *et al.* (2002). In this earlier investigation, reflection and transmission matrices comprised partitions corresponding to shear and dilatational waves and their cross coupling, while here the division is into Helmholtz and modified Helmholtz waves and their coupling.

For a particular grating, we define the matrices **R**^{±} and **T**^{±}, characterizing reflection and transmission of an incident plane wave coming from above or below the grating, respectively. It is necessary to specify the reflection and transmission coefficients *R*_{p} and *T*_{p} for Helmholtz type plane waves and and for the modified Helmholtz type waves. This must be done for the whole range of incidence channels enumerated by an integer *q*.

If *N* stands for the order of truncation of the multipole expansions, then the matrices **R**^{±} and **T**^{±} consist of four (2*N*+1)×(2*N*+1) blocks arranged as follows:
10.1
Here, the (*p*,*q*)th element of stands for the reflection coefficient *R*_{p} when the incident wave comes from above the grating, and its only non-zero amplitude component is for channel *q*. The (*p*,*q*)th element of contains for the same incident wave data. Correspondingly, the (*p*,*q*)th element of contains *R*_{p} when the incident wave comes from above the grating, and its only non-zero amplitude component is for channel *q*. Finally, the (*p*,*q*)th element of stands for the reflection coefficient with the same incident wave as for . The matrix **R**^{−} is filled in the same way, but with the incident wave now coming from below the grating, and the corresponding definitions are valid for the four blocks of the transmission matrices **T**^{±}.

Next, the reflection and transmission matrices we have introduced have a phase origin at the centre of the grating. When we use the grating as part of a stack, it is convenient to have such matrices with the phase origin half way between the successive elements. This is done using a diagonal propagation matrix . Let η be the distance between the gratings. If the phase origins for reflected and transmitted fields are placed at *y*=±(1/2)η, respectively, then the entries of the propagation matrix are
10.2
where if *u* corresponds to a plane wave of Helmholtz type and if *u* corresponds to a plane wave of modified Helmholtz type. The reflection and transmission matrices with the phase origins at *y*=±(1/2)η are then
10.3
which can be used for reflection and transmission matrices corresponding to incidence from above or below the grating.

Referring to fig. 3 of Platts *et al.* (2002), we suppose we know two matrices and characterizing reflection and transmission by a stack of *s* grating layers. We add a single layer characterized by reflection and transmission matrices and on top of the stack to give a new stack of *s*+1 layers. The reflection and transmission matrices of this new stack are given by the recurrence relations
10.4
and
10.5
The recurrence procedure is to fill and using the formulations of §§6 or 8 and deduce the reflection and transmission matrices for incidence from below ( and ) using the symmetry of the grating. These matrices are converted to phase-shifted form using equation (10.3) and used to start the recurrence. After *N*_{g}−1 applications of the relations (10.4) and (10.5), we have formed the desired matrices and characterizing the stack composed of *N*_{g} gratings. These are rephased to have their phase origins at the centres of the top and bottom layers of the stack using
10.6

## 11. Modes trapped in a pair of gratings

We now consider the reflection and transmission of a pair of identical gratings, each composed of clamped inclusions of zero radius. Our treatment generalizes that of Evans & Porter (2007), in that we include the effects of multiple plane waves of both Helmholtz and modified Helmholtz types in the interaction between the gratings, rather than just those of zero order.

Using equations (8.6)–(8.9), we can fill in the reflection and transmission matrices of the identical gratings, and thus the matrices padded to take into account the phase changes on propagation across half the distance η between the gratings: 11.1 The reflection and transmission matrices of the pair of gratings are then given by equations (10.4) and (10.5): 11.2

We show in figure 5 the energy transmitted as a function of *β* through a pair of gratings with unit separation for an angle of incidence *θ*_{i}=30^{°} (compare with figure 2 for a single grating). The Rayleigh anomaly (McPhedran & Waterworth 1972) associated with the onset of the order *p*=1 is now a very sharp spike in transmittance, followed by a relatively broad transmission resonance. Of course, the most intriguing feature of figure 5 is the very sharp transmission resonance at *β*=*β*_{t}=3.58221, where the transmittance attains 100%. The neighbouring half-power points occur at *β*=3.58187 and *β*=3.58253, giving a quality factor for the resonance of around 5400. This is extraordinarily large, given that the structure consists of only two gratings, but may be understood with reference to figure 2. Indeed, the resonant value of *β* corresponds in figure 2 to a point of reflectance quite close to 1 (in fact, 0.99955). Energy entering the region between the gratings is trapped owing to the high grating reflectance and adds up owing to constructive interference, with the plane wave amplitudes between the gratings exceeding the amplitudes outside the gratings by a large factor. The matrices giving plane wave amplitudes for the downward and upward directions between the gratings are in fact given by
11.3

Note that the determinant of the matrix attains its minimum value (around 0.0015) for *β* corresponding exactly to the value (3.58221) for which the transmittance of the grating pair is unity. The fact that this determinant value is small but non-zero shows that the value of *β* for which a true resonant solution occurs is in fact complex, with real part close to 3.58221, but with a small imaginary part.

As *β* increases above *β*_{t}, the transmittance falls rapidly and in fact has a zero at *β*_{0}=3.59194. This sort of behaviour of the reflectance or transmittance passing through the extreme values of zero and unity in a narrow frequency range is typical of Fano resonances (Fano & Cooper 1968), where a system has a rapid variation with frequency superimposed on a slower variation. In this case, the slower variation is due to the individual grating, as shown in figures 2, 4 and 5, while the rapid variation is due to the phase term incorporated in the propagation matrix of equation (10.2).

As the angle of incidence is changed from *θ*_{i}=30^{°}, the coincidence between the values of *β* at which the single grating reflectance is unity and the pair of gratings has the transmission resonance becomes less striking. For example, for *θ*_{i}=27^{°}, the energy transmitted for the grating pair and the single grating is shown in figure 6. The energy transmitted for the grating pair resonates at *β*=3.59999, while the maximum of the single grating reflectance now occurs at *β*=3.696. This less perfect alignment of these values results in a somewhat smaller quality factor of the resonance (now around 1700).

## Acknowledgements

R.C.M. acknowledges support from the Australian Research Council under its Discovery Grants Program and partial travel support from the Research Centre in Mathematics and Modelling of the University of Liverpool. C.G.P. acknowledges support from a start-up grant from the University of Technology, Sydney.

## Footnotes

- Received June 5, 2009.
- Accepted July 16, 2009.

- © 2009 The Royal Society