## Abstract

We study the scattering of elastic waves by platonic clusters in the time domain, both for plane wave excitations and for a specified initial wave profile. We show that we can use an analytical extension of our problem to calculate scattering frequencies of the solution. These allow us to calculate approximate solutions that give the flexural wave profile accurately in and around the cluster for large times. We also discuss the early-time behaviour of flexural waves in terms of the classical models of Sommerfeld and Brillouin.

## 1. Introduction

We present a solution in the time domain to the problem of a vibrating elastic thin plate, which is pinned at a finite number of points. These points form a finite subset of a regular lattice or crystal, and we will use the term *cluster* to describe such a system. The description *platonic cluster* then refers to a vibration problem for a thin plate governed by the biharmonic equation, with the plate having constraints on its motions, or inclusions imposed at a cluster of points.

The problem we discuss was studied in the frequency domain, initially by Evans & Porter (2007) and recently by McPhedran *et al.* (2009) and Movchan *et al.* (2009). McPhedran *et al.* (2009) highlighted the simply analytical form of the vibration band structure associated with arrays of pinned points and the form of defect modes associated with a finite set of points. We consider here the time-domain solution constructed from the frequency-domain solution. The time-domain solution gives new insights into the physics of the system, especially those connected with the presence of near-trapping behaviour. The solution method is based on a generalized eigenfunction expansion, in which the solution is expanded in terms of the frequency-domain solutions for incident plane waves. This method could be extended to more complex scattering problems than just pinned points, but our problem has the advantage that it is especially simple to solve in the frequency domain, and the physical behaviour of the elastic wave interactions emerges clearly.

What we will demonstrate, using simulations and by analytical methods, is that the elastic wave interaction with various examples of clusters divides into three time regions. In the first, multiple wave interactions with the structure have not had sufficient time to develop, and the front of the wave travels with the same velocity as in a uniform plate. In the second interval, very complicated wave patterns arise, as multiple resonant wave patterns compete with each other. In the last interval, the longest-lived resonant mode dominates, and its wave pattern gives that for the cluster. It is in this last interval that the approximate solution we develop becomes most accurate and useful.

The physics we demonstrate using the case of elastic waves in thin pinned plates is generic to a wide class of temporal problems involving waves in structured media. It is advantageous to study the problem discussed here since it gives rise to interesting and clear cut wave-trapping behaviour, evident in even small clusters. It is also very simple to implement numerically, and a wide range of geometries and wave phenomena in addition to those studied can easily be investigated using the same techniques developed here.

The generalized eigenfunction method is an extension of the standard eigenfunction method for compact self-adjoint operators to non-compact self-adjoint operators. The method was developed originally for Schrodinger's equation by Povzner (1953) and Ikebe (1960) and has since been extended to various other equations, including the wave equation, as discussed in Wilcox (1975). The method is quite general and can be applied to many different problems, for example, it has been applied extensively to water waves (Hazard & Lenoir 2002; Hazard & Loret 2007*a*; Hazard & Meylan 2007; Meylan 2009; Meylan & Eatock Taylor 2009). However, the method has only recently been developed as a numerical method (Meylan 2002; Hazard & Meylan 2007; Meylan 2009; Meylan & Eatock Taylor 2009; Peter & Meylan 2010; Fitzgerald & Meylan 2011).

The singularity expansion method (SEM) is a method to approximate the time-dependent equations by an expansion using the scattering frequencies or resonances. The scattering frequencies are the singularities of the analytical continuation of the resolvent. The SEM is usually derived from the Laplace transform solution (Baum 1976), however, a method to connect the SEM with the generalized eigenfunction method was found by Meylan & Eatock Taylor (2009) in the context of the two-dimensional wave equation (which was then applied to model water-wave scattering). This new SEM has several advantages over the traditional method based on the Laplace transform, and it gives the expansion formula directly.

The outline of the paper is as follows. In §2, we derive the governing equations in the time domain. In §3, we develop the equations in the frequency domain. In §4, we show that the solution in the frequency domain can be approximated from knowledge of the scattering frequencies. In §5, we show how we can solve the equations in the time domain for an incident wave packet and for an arbitrary initial condition. In both cases, we show how we can approximate the solution using the scattering frequencies. In §6, we present numerical results, with both figures and animations being used to illustrate the evolution in time of waves and pulses interacting with clusters of pinned points of various sizes. We finish with a brief summary in §7, and an appendix giving comments on the wave phenomena evident in each of the animations.

## 2. Governing equations in the time domain

We begin with the equation for an elastic plate in the time domain
2.1where *τ*(**x**,*t*) is the displacement of the plate, *t* is the time and **x**=(*x*,*y*) represents a point in the plate in terms of Cartesian coordinates in two dimensions. Note that we have assumed that the variables are non-dimensionalized so that the ratio of the mass and stiffness terms is unity (by scaling either the time or length). We consider here the simplest scattering problem: that of a finite number of pinned points at which the displacement is zero. This is a simple problem to solve in the frequency domain and it has been the subject of considerable recent work, owing to the interesting behaviour of its solution, especially in the context of periodic arrays (Evans & Porter 2007; Movchan *et al.* 2009). Consequently, we specify the displacement at a set of *N* points **x**_{j} to satisfy the condition
2.2The solution of equation (2.1) with boundary conditions given by equation (2.2) requires conditions in the time domain. We consider two types of conditions. The first is that for times with a large negative part, we have an incident plane wave packet and this condition is presented in §5. The second is to prescribe initial conditions at *t*=0 of the form:
2.3aand
2.3bWe note that the theory which we develop here could be extended to more complex scattering problems, for example, circular inclusions as treated in McPhedran *et al.* (2009).

## 3. Equations in the frequency domain

We present here the equations in the frequency domain and show how they can be easily solved by using the fundamental Green function. The frequency-domain equations are found by assuming time harmonic dependence of the form e^{−iωt}, i.e.
3.1Consequently, equation (2.1) in the frequency domain becomes
3.2and the boundary conditions are unchanged, i.e.
3.3To solve the frequency-domain equation, rather than imposing initial conditions, we need to consider the solution at infinity. We impose the condition that the solution consists of an *incident wave* and a *scattered wave*. The equation in the frequency domain can be factorized into a Helmholtz and anti-Helmholtz (or modified Helmholtz) equation. For this reason, the scattered wave satisfies the standard Sommerfeld radiation condition for the Helmholtz equation plus a condition that the solution is bounded at infinity. We consider the system excited by an incident or imposed wave, which we denote by *η*^{I}. In what follows we consider two types of incident waves. The first is a plane incident wave of the form:
3.4where is the wavenumber and *χ* is the incident angle of the plane wave (with respect to the Cartesian axes and measured from *Ox* in the standard way). These waves are used to simulate an incident wave packet. We also consider an imposed cylindrical wave of the form:
3.5where *J*_{n} represents a Bessel function of the first kind (Abramowitz & Stegun 1970) and (*r*,*θ*) are polar coordinates. The cylindrical waves will be used later in our generalized eigenfunction expansion.

The solution to the scattering problem of a set of pins is given by writing the solution as a combination of incident waves and Green functions as done by Evans & Porter (2007). We define the total displacement of each incident wave to be
3.6or
3.7with
3.8where is a Hankel function of the first kind (Abramowitz & Stegun 1970) and *g*_{k} is the fundamental Green function for this problem, which satisfies
3.9where *δ* represents the Dirac delta function. To find the coefficients *a*_{j} (which will be different for different incident waves or pin positions), we simply solve the problem
3.10which is an *N*-dimensional matrix equation. Note that *η*^{I} is used to mean either or . This equation can be written as a matrix equation for the unknown vector **a** of coefficients *a*_{j} as
3.11where **G** is the matrix that follows from the left-hand side of equation (3.10) and **f** is the vector that follows from the right-hand side. In what follows we will use the notation **f**_{χ} and **f**_{n} to distinguish the two types of incident waves.

## 4. Scattering frequencies

Equation (3.11) is the equation to calculate the scattering for real frequencies. The solution to this is given by
4.1We now consider equation (3.11) for complex *k*. We find numerically that the matrix **G**(*k*) has zero eigenvalues for a set of points *k*=*k*_{p} in the lower half of the complex plane. The fact that they only occur in the lower half complex plane is exactly what is expected theoretically from causality, and from our choice of e^{−iωt} time dependence (with positive *ω* corresponding to positive *k*). From equation (4.1) for *k*=*k*_{p}, the solution **a** will have a pole, i.e. the system will have infinite response. The values *k*_{p} are called scattering frequencies and they play a role similar to eigenvalues. We will show that they can be used to calculate the solution in the time domain for large times.

In general, there is no requirement for the analytical continuation to be meromorphic. However, for the elastic plate equation, by analogy with the wave equation, we would expect the analytical continuation to be meromorphic, excepting the branch cut, which arises from the logarithm in the fundamental solution (i.e. from the Hankel function). This is exactly what we find numerically. We will place the branch cut on the imaginary axis. It is possible to approximate the solution in the time-domain by deforming the contour of integration in the complex plane and including the contribution from the singularities explicitly. This method is known as the SEM. We will use a modification of the SEM in what follows.

For the purpose of our numerical results, we consider three cases. The first is 4 pins arranged in a square with spacing 2. The second is a 3 by 3 cluster of 9 pins, also with spacing 2 and the third is a 4 by 4 cluster of 16 pins with spacing 2. In all cases, the pins lie on a grid and all the clusters are centred at (0,0). We refer to these as the 4, 9 and 16 pins cases, respectively. Figure 1 shows the calculation of the roots using a complex search algorithm (Meylan & Gross 2003) with the single roots denoted by a cross and the double roots denoted by a circle (note that these arise owing to the symmetry of our problem). For the 4-pin case, there are two complex scattering frequencies close to the real axis. These separate into multiple scattering frequencies for the 9 and 16-pins case. In what follows, we focus on the scattering frequencies with real part close to *k*_{p}=1.8. The single value for 4 pins breaks into four scattering frequencies for the 9-pin case and nine for the 16-pin case (counting double roots twice).

### (a) Calculation of the residues

Suppose we have a scattering frequency at a complex wavenumber *k*_{p}. The scattering frequency *k*_{p} is associated with an eigenvector **u**_{kp} with the property that
4.2Near a simple pole, we can write the inverse of **G**(*k*) as
4.3where **A**(*k*_{p}) is the residue which is connected with a projection onto the eigenspace associated with **u**_{kp} and **v** is an arbitrary vector. The expression for **A**(*k*) follows from Hazard & Loret (2007*b*, theorem 3.5, p. 934) using a result owing to (Steinberg 1968) when the eigenspace has dimension one. The case of a double root will be presented separately for reasons of clarity. In practice, the double roots situation is important because double roots do occur in problems with symmetry (which we consider here). It can be shown that near the simple root *k*_{p},
4.4where is the eigenvector of the adjoint of **G** with eigenvalue zero (which occurs at *k*^{*}_{p}). The dot denotes vector product and the star denotes complex conjugate. The action of the right-hand side is to multiply a column vector by the row vector to produce a scalar and then to multiply this by the column vector **u**_{kp} normalized by the denominator (which is also a scalar). **G**^{(1)} is the derivative of **G** at *k*_{p} given by
4.5

#### (i) Double root formula

In the case of double roots, we assume that there are two linearly independent eigenvectors associated with *k*_{p}. The case where there is only a single eigenvector would definitely be exceptional, and we have not encountered this situation in any numerical results. Such a situation would also imply more complicated time-dependent behaviour, corresponding to a non-physical problem.

We let **u**_{kp} and **v**_{kp} be two linearly independent orthogonal eigenvectors, and we choose two eigenvectors of the adjoint, which are given by and with the properties that and the added requirement that they are orthogonal. Then
4.6As mentioned previously, we are assuming the double root corresponds to two distinct eigenvectors (i.e. there is not a Jordan block). We do not know of a mathematical explanation of why this has to be the case, but the behaviour in the time domain would include a secular term if we have a Jordan block, which seems highly unphysical.

#### (ii) Mode shape

We can find the mode shape (which is similar to an eigenvector) associated with each of the scattering frequencies from
4.7where *u*_{j} is the *j*th component of **u**_{kp}. Note that the function *U*_{kp}(**x**) is not bounded as because the Hankel function of the first kind grows for complex argument with large negative imaginary part. It might seem that using a function that grows to approximate a solution, which decays, is contradictory. However, the approximation is only valid for large times and small distances. While the Hankel function grows with distance, there is also an exponential decay in time. It follows that near the simple scattering frequency
4.8This formula can be extended to the case of double roots straightforwardly using the two parts of equation (4.6).

Figure 2 shows the mode shape associated with each of the two resonances close to the real axis for the 4-pin cluster. The value for the complex scattering frequency is given for each figure. Figure 3 shows the four mode shapes associated with the scattering frequencies close to 1.80 for the 9-pin cluster. Note that the single mode for 4-pins has split into four modes for the 9-pin cluster. Note also that the double root gives rise to two mode shapes for the same frequency value and that these mode shapes are related by a simple rotation. For the 16-pin cluster, the single mode for 4-pins has now split into nine modes as shown in figure 4. The splitting of the mode follows from group theory, as described in McIsaac (1975).

## 5. Solution in the time domain

We show here how the time-dependent equations can be solved using the single frequency solutions for two kinds of initial conditions. The first is an incident wave packet from infinity, and the second is for prescribed arbitrary initial conditions. In both cases, we present approximate solutions using the scattering frequencies.

### (a) Incident plane wave packet

The frequency-domain solutions can be combined to give solutions in the time domain, although this problem is more subtle than is often believed. The simplest problem is to consider a wave packet, which is composed of plane waves travelling in the same direction, but with different frequencies. In the absence of scattering, this problem is trivial, and we write
5.1where we have used the relationship *ω*=*k*^{2}. We assume that packet at *t*=0 is given by a function , with *χ* specifying the propagation direction of the packet. Therefore
5.2and can be calculated from the Fourier transform formula
5.3If we now consider the case when there is scattering, and excite the system with the same incident wave packet as , we can calculate in exactly the same way using equation (5.3), and we write
5.4

### (b) Approximation using scattering frequencies

We can use our results for scattering frequencies to approximate the solution for long time. The formula which we derive relies on a deformation of the contour of integration. We then assume that all terms, besides the contribution from the calculated scattering frequencies, are small and can be neglected. This assumption works well in the numerical examples where the finite cluster has a small number of scattering frequencies very close to the real axis, and where we excite the system close to the frequency of the scattering frequencies. We do not claim that it works well in all situations and we do not have formulae to give an expression for the size of the terms we have neglected. However, for situations such as we consider here, where there is near trapping of energy these expansions work very well.

The path of integration that we consider is given by the solid line along the real axis shown in figure 5. This integral can be evaluated instead by deforming the contour of integration and determining the integral shown by the dotted line closed in the lower complex plane as well as the contributions from the singularities denoted by the circles. Our approximation is based on neglecting all the terms denoted by dashed lines, i.e. we do not include contributions from the integral path in the lower complex plane nor from singularities far from the real axis.

We proceed as follows. We substitute the approximation (4.8) into equation (5.4), allowing for *P* scattering frequencies (counting double roots twice), and excluding the term owing to the plane incident wave. The plane incident wave part of the solution has a particularly simple form and we assume that we are interested only in the solution after it has passed by our finite cluster. This gives us
5.5where **f**_{χ} denotes the right-hand side of equation (3.11) given by a plane incident wave. We can then close the contour in the lower half plane. This gives rise to contributions from the scattering frequencies and also from the contours of integration in the lower complex plane. We assume that the contribution from the contours of integration can be neglected. This assumption seems reasonable for the contour with large *k* because either will be small or the exponential part will vanish. However, there will be some contribution from the integral along the imaginary axis and this contribution will not vanish, and may include algebraically decaying parts. However, for our examples, is assumed to be small for small *k* (since the incident wave packet is localized) and subsequently the contribution from this integral can also be expected to be small. Under these assumptions, we obtain
5.6

### (c) Generalized eigenfunction expansion

The generalized eigenfunction expansion method is based on expanding the time-domain solution using the frequency-domain solution in a similar manner to before. However, in the previous section, we considered an incident plane wave packet, which considerably simplified the problem and allowed us to compute the solution in the time domain directly, using the standard Fourier transform (when going from the spatial to the frequency domain). We considered this case separately because it is of significant physical interest, and it corresponds to a standard scattering problem. However, if we consider an initial value problem, in which the value of *τ* and its time derivative are prescribed at *t*=0, then we need to use a generalized eigenfunction expansion in which we have to use a special transform to go from the spatial to the frequency domain. The generalized eigenfunction expansion requires a complete set of imposed waves (i.e. we cannot consider waves which are incident from a single direction) and we calculate it using the cylindrical waves *η*_{n} and summing over *n*. This choice was also made by Meylan & Eatock Taylor (2009). An equivalent formula using plane incident waves from all directions could be developed. We can write the solution in the time domain as a superposition of these waves, i.e.
5.7It is clear that this is a solution of equation (2.1), however, we still need to find the values of and so that at *t*=0 we can satisfy the initial conditions (2.3). This requires that
5.8aand
5.8bThe generalized eigenfunction method rests on one key idea, that the inner product of the eigenfunctions normalizes the same way with and without the scatterers. The inner product needs to be chosen so that the operator (here Δ^{2}) is self-adjoint. For our case, this is just the standard inner product
5.9where the star denotes the complex conjugate. Using the results that the eigenfunctions have the same normalization with and without the scatterers (in which case *η*_{n}=*η*^{I}_{n}), we obtain
5.10where *δ*_{nn′} is the Kronecker delta. Therefore, taking inner products with respect to *η*_{n}(*r*,*k*) in equations (5.8a) and (5.8b), we obtain
5.11aand
5.11b

### (d) Approximation using scattering frequencies

We can approximate the solution in the time domain using the approximation (4.8) and closing the contour in the complex plane, again ignoring all other contributions except those owing to the scattering frequencies. The justification for this is exactly as before. We also assume *h*(**x**)=0 for simplicity and that any cross contribution between the scattering frequencies is negligible (from their near orthogonality). This gives us
5.12where **f**_{n} denotes the right-hand side of equation (3.11) given by a cylindrical wave of order *n*.

## 6. Results

We now present some results, beginning with the response of our clusters to a plane incident wave packet. The brief comments in this section highlight some principal features of the five animations accompanying this paper, with a more extended description being given in the appendix. The animations themselves are central to understanding the interaction of waves with the cluster, and we recommend that the reader carefully examine these, since they contain more information than can be compactly described here.

We consider a Gaussian plane wave packet centred at frequency *k*_{a} of the form:
6.1Since the medium is highly dispersive, this corresponds to an incident packet which peaks at *t*=0 and *x*=0 (in the absence of scattering) exactly at the centre of the cluster with amplitude 0.5. We fix the incident direction to be *χ*=0, so that the wave propagates in the positive *x*-direction. Figure 6 shows the response of the cluster for a normally incident wave packet with *k*_{a}=1.79. Shown in the figure is the displacement *η* for times *t*=−0.5, 0.5, 1.5 and 2.5 as a three-dimensional plot, and also as 2 line plots along the lines *y*=0 and *y*=2. The *exact* solution is the solution calculated using equation (5.4). While this solution is calculated numerically, the solution is exact up to the accuracy of the numerical method and this is the reason we label this solution as exact. Also shown in the line plots are the *approximate* solutions calculated using equation (5.6). This solution is clearly approximate and will only be valid for large times, as can be seen in the results. We also show a third curve, the *undisturbed* wave, which is the wave profile that would be present if the scatterers were absent. The approximate solution is based on the scattering frequencies close to 1.79. Figures 7 and 8 show the equivalent plots for 9 and 16 pins. These figures are snapshots from the electronic supplementary material, movie files, which show the solution in greater detail.

The initial physics of the interaction is made clear in figures 6–8. The front of the pulse arrives simultaneously at points with the same value of *x* (the direction of propagation of the plane wave pulse), regardless of whether they are inside or outside the cluster. This is consistent with the classical analysis of Sommerfeld and Brillouin, which is well described by Stratton (1941) and Jackson (1975). To quote Stratton (p. 337), with *c* denoting the wave speed in free space,
‘the wavefront velocity is always equal to

*c*, no matter what the medium.’

Obviously our situation is more complicated since the medium is highly dispersive, but our pulse is relatively narrow banded. To summarize, the wave interactions with the pins take a definite time interval to develop and stabilize and during this time interval the interaction of the wave with the cluster is weak.

The wavefront is also referred to as the *Sommerfeld precursor wave*. Note that a second precursor wave is associated with the name of Brillouin, and it occurs between the Sommerfeld wave and the achievement of a steady-state response. We have investigated that this *Brillouin precursor wave* is not evident in figures 6–8 and also in other examples, and its absence may be explained by the fact that Brillouin's analysis relied critically on the wave properties of the medium being represented accurately by a dispersion relation with a single pole and zero (together with their negative frequency counterparts). This model clearly cannot represent our wave bandgap system, as is evident from figure 1. Furthermore, the highly dispersive nature of our medium can be expected to smooth out the response.

The behaviours at times between the arrival of the Sommerfeld precursor and the onset of the asymptotic behaviour are quite variable with position, and with the pulse properties. A cluster of moderate-to-large size can support a number of modes. These will have shorter lifetimes if they have antinodes close to the edge of the cluster, and longer lifetimes if antinodes are deep in the cluster. The former are quickly excited, but decay rapidly. The latter take more pulse cycles to be excited, but can have much larger amplitudes and longer decay times in consequence. The competition between these modes gives rise to a complicated range of temporal and spatial dependencies for intermediate times which cannot be given any simple mathematical description, beyond the generalized eigenfunction formula equation (5.7). However, once the incident pulse has passed through the cluster, the approximation (5.12) quickly gives an accurate representation of the vibration patterns in each of the first three animations. We note that the accuracy of this approximation does not require that the system has a single mode character. This means that the corresponding pattern of vibration can be somewhat complicated for the case when we have a large number of pinned points.

We now consider the system excited by an initial displacement. This case is computationally more challenging. The initial displacement was chosen to be of the form: 6.2

Figure 9 shows the displacement for an initial pulse with *x*_{p}=*y*_{p}=1 for the 9-pin cluster and figure 10 is the equivalent plot for the 16-pin cluster except that *x*_{p}=*y*_{p}=2. This time we use equation (5.7) to calculate the exact solution and equation (5.12) to calculate the approximate solution. Again, we note that by the exact solution, we mean a numerical evaluation of equation (5.7). It is interesting to note that we can obtain information about the numerical accuracy of our solution by how well it approximates the imposed initial condition at *t*=0. As before, the short-time behaviour is complicated and cannot be given a simple approximation. However, once the wave has sufficient time to propagate through the cluster, we get an accurate approximation using the small number of modes as before. Again, we have a competition between resonant modes so that only in the fourth animation does a stable vibration pattern develop. In the more complicated situation shown in the fifth animation, we see a pattern formed by the superposition of a number of modes.

## 7. Summary

We have presented a solution to the problem of wave scattering from pinned elastic plates in the time domain. The simple nature of the pinned boundary conditions means that we can solve the problem very efficiently in the frequency domain, and we use the frequency-domain solution to construct the solution in the time domain. We consider two cases, the first in which the system is excited by an incident plane wave pulse and the second in which the system is excited by an initial displacement of the plate. We have also presented an approximate solution based on the complex scattering frequencies. This simple approximation has proved to be very accurate for large times. Numerical results showed that, after an initial transient period, the solution was well represented by the approximate solution.

## Acknowledgements

M.H.M. acknowledges the support of the University of Auckland Staff Research Grants. R.McP. acknowledges the support of the Australian Research Council through its Discovery Grants Scheme. The authors would like to thank Dr Garry Tee and Mr M. Smith for their careful reading of the manuscript as well as the reviewers for their helpful feedback.

## Appendix A

#### (a) Description of the animations

##### (i) Animation

In this simulation, a Gaussian wave packet approaches a cluster of 4 pins, travelling parallel to the *x*-axis. The displacement *τ* is given as a function of position for increasing times. In the two panels on the right, line plots of the displacement as a function of *x* are given for various times, for two values of *y*: *y* = 0 in the centre of the cluster, and *y*=2 outside the cluster. Three lines are given in each panel: the exact (numerical) solution, the approximate solution (based on a single mode in equation (5.12)) and the undisturbed solution (where there are no pins). It is clear that the approximate solution varies in an exponential fashion with time, and only begins to be accurate when the wave packet has entered the cluster.

Note that for time *t*<0, the front of the exact solution and that of the undisturbed wave travel together (with both exiting right for *t*≈0.5). The centre of the pulse reaches the cluster at *t*=0, and begins to be strongly reshaped. At *t*≈1.4, the approximate model begins to resemble closely the exact solution around *x*=0. At this time, the wave packet has largely left the centre of the cluster. The region of accuracy for the approximate solution then grows with increasing time, as the wave is increasingly driven by the dominant mode. By *t*=3.6, the approximate solution works well in the whole computational domain. The dominant mode slowly decays in amplitude, but with constant shape as *t* increases. Note that its shape for these values of *t* along *y*=2 (outside the cluster) replicates that at *y*=0 (inside the cluster).

##### (ii) Animation

In this animation, the same Gaussian pulse interacts with a cluster of 9 pins. The approximate solution is now found using four modes (although only two of these are excited by waves incident normally). As before the exact and undisturbed wave fronts travel together and leave the computational domain at *t*≈−0.5. We also have a period of complex scattering as the pulse interacts with the 9 pins. The approximate solution does not become accurate at the centre until *t*≈1.3 (note that it is still not accurate for the line *y*=3 at this time). At around *t*≈3.5, the approximate solution is accurate for the entire computational domain. These times are the same as for the 4-pin animation because the wave packet peak reaches the centre of the cluster at *t*=0 for both cases. As mentioned previously, only two modes are excited by this particular incident wave. The two modes oscillate with different frequencies and they decay at different rates. This leads to a complex beating effect, so that energy moves from the front to the back of the cluster and vice versa. For very long times, the mode with the slowest decay will dominate, but we do not show the solution for the very large times where this would become apparent.

##### (iii) Animation

In this animation, the same Gaussian pulse interacts with a cluster of 16 pins. The results are again very similar to those for the 4- and 9-pin clusters. In this case, nine modes are used in the approximate solution, although only five are excited by the normally incident wave.

##### (iv) Animation

Animation 4 shows the solution for 9 pins for an initial Gaussian at *t*=0 centred at (1,1). In this figure, we no longer have an undisturbed wave and the first part of the solution for negative *t* may be thought of as showing the initial wave, which would be incident to excite this Gaussian initial condition. The solution in the time domain is symmetric about *t*=0 while the approximate solution, which is calculated from four modes with frequency around *k*=1.8, does not have this symmetry. The initial condition has been chosen to excite these modes. The approximate solution only becomes valid after *t*≈1 at the point (1,1) and quickly becomes valid by *t*≈2 throughout the computational domain. The high-frequency waves that are excited by the initial condition can be seen quickly dissipating as they flow away from the cluster.

##### (v) Animation

Animation 5 shows the same Gaussian pulse, this time centred at (2,2) for the 16-pin array. In this case, the evolution of the vibration pattern in the cluster is more complicated. The displacement peak will evolve to be further from the edge of the cluster, so that it will persist for large times (since strong vibrations near the edge lead to rapid decay). Such a peak movement requires excitation of many of the nine modes, and during the time interval considered the peak actually *wanders* around the cluster. The vibration pattern has not stabilized by *t*=10, however, the approximation (5.12) is still an accurate representation of the numerical solution.

- Received April 14, 2011.
- Accepted June 29, 2011.

- This journal is © 2011 The Royal Society