## Abstract

We consider the propagation of free surface waves on an elastic half-space that has a localized geometric inhomogeneity perpendicular to the direction of wave propagation (such waves are known as topography-guided surface waves). Our aim is to investigate how such a weak inhomogeneity modifies the surface-wave speed slightly. We first recover previously known results for isotropic materials and then present additional results for a generally anisotropic elastic half-space assuming only one plane of material symmetry. It is shown that a topography-guided surface wave in the present context may or may not propagate depending on a number of factors. In particular, they cannot propagate if the original two-dimensional surface wave on a flat half-space is supersonic with respect to the speed of anti-plane shear waves. For the case when a topography-guided surface wave may exist, the existence and computation of wave speed correction is reduced to the solution of a simple eigenvalue problem whose properties are previously well understood. As a by-product of our analysis, we deduce that there exists at least one topography-guided surface wave on an isotropic elastic half-space, and that it is unique when the geometric inhomogeneity has sufficiently small amplitude.

## 1. Introduction

This paper is concerned with the effects of a localized hump or depression on the propagation of an elastic surface wave along the otherwise flat surface of an elastic half-space; see figure 1. A (classical) surface wave is a wave propagating along a free flat surface, for which the associated displacement and stress decay away from the surface. They appear in seismology, signal processing, non-destructive testing and within many other physical scenarios. Lord Rayleigh [1] was the first to show that the plane traction-free surface of an elastic isotropic half-space could support such a surface wave. The intuitive approach initially adopted for an isotropic half-space is not capable of resolving the existence and uniqueness question for surface waves on a generally anisotropic elastic half-space. For the latter, the Stroh formalism [2] is the necessary tool and it was with this tool that Barnett *et al.* [3], Barnett & Lothe [4] and Chadwick & Smith [5] gave the first existence and uniqueness proof. We refer to Barnett & Lothe [6] for an updated account of the existence and uniqueness theory and to the book by Ting [7] for a more detailed description of the theory and for a comprehensive collection of relevant references.

A surface-wave solution may also be referred to as a trapped mode in the sense that propagation energy is confined to within a few wavelengths of the free surface. At the same time as extensions of Rayleigh's classical surface-wave solution were made to anisotropic elastic and/or viscoelastic materials, much interest was also shown in trapped modes around thin rectangular plates, wedges, or similar structures attached to a half-space. There are two important limits corresponding to this set-up. When the rectangular plates are thin and long, the trapped modes are expected to localize near the free edge and the associated modes are also known as edge waves; e.g. Konenkov [8], Norris *et al.* [9], Fu & Brookes [10], Lawrie & Kaplunov [11], and references therein. When the extrusions are almost flat, we expect the trapped modes to be the classical surface-wave mode slightly modified by the geometrical inhomogeneity. It is the latter that we are interested in the present paper. For the intermediate case, existence results have been obtained by Bonnet-Ben Dhia *et al.* [12] using the min–max principle, and numerical results have been given by Burridge & Sabina [13]; see also Biryukov *et al.* [14] for a review of some approximate methods.

Our present work is partly motivated by the recent series of studies by Kaplunov *et al.* [15], Gridin *et al.* [16,17], Adams *et al.* [18] and Postnova & Craster [19,20], in which a multiple-scale approach was fruitfully used to study trapped modes in a variety of problems, but attention was invariably focused on isotropic materials. In this paper, we demonstrate that the same methodology, coupled with the Stroh formalism, can be used to deal with anisotropic materials. It turns out that our general formulation also enables us to settle the existence and uniqueness question left open by Adams *et al.* [18] concerning topography-guided surface waves on an isotropic elastic half-space.

The paper is organized as follows. After formulating the problem in §2, we present and solve the leading-order and second-order problems in the subsequent two sections, respectively. In §5, we write down the third-order problem and obtain the amplitude equation by imposing a solvability condition. It is shown that this amplitude equation can be reduced to the time-independent Schrödinger equation whose spectral properties are well known. After demonstrating in §6 that our formulation can recover the isotropic results of Adams *et al.* [18], we solve our amplitude equation in §7 numerically. Illustrative examples are given that demonstrate a rich variety of behaviour associated with anisotropy: we may have zero, a single or multiple topography-guided solutions, and the speed may be higher or lower than the speed corresponding to a flat surface. The paper is concluded in §8 with additional discussions.

## 2. Problem formulation

We consider surface-wave propagation on a half-space whose surface is defined by
2.1where *x*_{i} are Cartesian coordinates, *h* is a smooth localized function and *ε* is a small positive parameter. We shall use the wavelength *L* as the reference length scale and assume that *h*/*L* is of order unity. The specified function *h* is thus slowly varying and is localized in the *x*_{3}-direction; see figure 1. The surface-wave problem is governed by the equation of motion and constitutive relation
2.2together with the traction-free boundary condition and decay condition
2.3where *u*_{i} and *σ*_{ij} (*i*,*j*=1,2,3) are the displacement and stress components, *ρ* is the density, *c*_{ijkl} are the elastic moduli and (*n*_{i}) is the unit normal to the free surface. Here and hereafter, we use the summation convention on repeated suffices and use a comma to signify partial differentiation with respect to spatial coordinates, and a superimposed dot to denote differentiation with respect to time.

We assume that the moduli obey the symmetry relations *c*_{ijkl}=*c*_{klij}=*c*_{jikl}, and the strong convexity condition
2.4For numerical illustrations, we shall consider two representative anisotropic elastic materials. The first is a transversely isotropic elastic material with moduli given by
2.5where *λ*,*α*,*β*,*μ*_{t} and *μ*_{l} are material constants and *m*_{i} is the preferred direction; see Spencer [21]. These constants are related to Young's moduli *E*_{l} and *E*_{t} and Poisson's ratio *ν*_{lt} as follows:
where . We shall take *ρ*=1852 kg m^{−3},*ν*_{lt}=0.324 and
2.6given by Rikards *et al.* [22] for a typical glass–epoxy composite. The second material is the cubic material silicon with material constants taken from Farnell [23],
2.7

For a surface wave propagating in the *x*_{1}-direction with speed *v*, the dependence on *x*_{1} and *t* is through *x*_{1}−*vt*, and so equations (2.2) and (2.3)_{1} may be written as
2.8and
2.9The free surface is defined by *x*_{2}+*h*(*εx*_{3})=0, and so its unit normal is in the direction (0,1,*εh*′(*εx*_{3})), where the prime indicates differentiation with respect to the argument *εx*_{3}. Thus, the traction-free boundary condition (2.9) may be written as
2.10We now introduce the variable transformation
so that in terms of the new coordinates, the free surface is given by *x*_{2}′=0. We have
and in terms of the new coordinates, (2.8) and (2.10) become
2.11and
2.12where here and hereafter, Greek subscripts range from 1 to 2, and a comma now signifies partial differentiation with respect to the primed coordinates. To simplify notation, we shall from now on drop the primes on the coordinates.

We now look for a perturbation solution of the form
2.13where *X*_{0} and *X*_{1} are constants and **u**^{(k)} (*k*=0,1,2,…) are functions to be determined. Obviously, the solution (*X*_{0},**u**^{(0)}) describes a surface wave when the surface is flat. Our aim is to find the leading-order correction *ε*^{2}*X*_{1} when the surface is specified by (2.1).

On substituting (2.13) into (2.11) and (2.12) and then equating coefficients of like powers of *ε*, we obtain, from the first three orders, respectively, that
2.14
2.15
2.16
2.17
2.18
and
2.19

## 3. The leading-order problem

The leading-order problem is the classical surface-wave problem, the solution of which may be written in the form
3.1where c.c. denotes the complex conjugate of the preceding term, *k* is the wavenumber and *f*(*kx*_{3}) is a function to be determined. The depth variation * z*(

*kx*

_{2}) may be represented by 3.2where

*E*is a 3×3 matrix whose eigenvalues determine the decay rate; see Fu & Mielke [24]. We shall initially set

*k*=1, and then show at a later stage (see the discussion leading to (5.16)) how to obtain the corresponding results when

*k*is not unity.

On substituting (3.1) with *k*=1 into (2.14) and (2.15), we find that the equations of motion reduce to
3.3and the boundary conditions to * t*(0)=0, where the reduced traction vector

*is given by 3.4and the matrices*

**t***T*,

*R*and

*Q*

^{(v)}are defined by their components 3.5On substituting a trial solution of the form

*=*

**z***e*

**a**^{ipx2}into (3.3), we obtain 3.6so that the values of

*p*are determined by setting the determinant of coefficient of

*to zero. For*

**a***v*in the subsonic interval [5], none of the values of

*p*can be purely real. Furthermore, when the six values of

*p*are distinct (as three pairs of complex conjugates), the surface-wave solution can be written as 3.7where

*p*

_{1},

*p*

_{2}and

*p*

_{3}are the eigenvalues of (3.6) with positive imaginary part,

**a**^{(k)}(

*k*=1,2,3) are the associated eigenvectors,

*=(*

**c***c*

_{1},

*c*

_{2},

*c*

_{3})

^{T}is a constant vector, The last expression in (3.7) then shows that in this case, the matrix

*E*in (3.2) is equal to −i

*A*〈

*p*〉

*A*

^{−1}.

When the six values of *p* are not distinct, the matrix *E* can be determined as follows. First, with the use of (3.3) and (3.4), the first-order derivatives * z*′(

*x*

_{2}) and

*′(*

**t***x*

_{2}) can easily be written as a linear combination of

*(*

**z***x*

_{2}) and

*(*

**t***x*

_{2}), thus leading to 3.8the Stroh formulation, with

*N*

_{1},

*N*

_{2}and

*N*

_{3}given by 3.9Next, we define the

*surface-impedance matrix*

*M*through 3.10Since the half-space is homogeneous, the above relation implies that 3.11On substituting (3.11) and (3.2) into (3.8), we find that the matrix

*E*in (3.2) and the surface impedance matrix are related by 3.12and that

*M*satisfies the Riccati matrix equation 3.13It follows from the definition (3.10) of the impedance matrix that the traction-free boundary condition

*(0)=0 can be satisfied only if 3.14which is the secular equation determining*

**t***X*

_{0}. When a root of this equation is found, a solution of

*(0) can be found by solving*

**z***M*(0)=0, and the corresponding surface-wave solution is then given by (3.1) and (3.2).

**z**Finally, we note that although (3.13) has multiple solutions, it has a unique solution for *M* that is positive definite for *v* less than the unique surface-wave speed, and this unique solution is what should be selected in calculating the surface-wave speed according to (3.14) and *E* according to (3.12) [24].

For an orthotropic elastic half-space whose axes of symmetry coincide with the coordinate axes, the three matrices *T*,*R* and *Q*^{(v)} take the simple form
3.15It can be shown that the impedance matrix *M* then assumes the form
3.16The Riccati matrix equation (3.14) can be solved exactly to find the unique solution mentioned above,
3.17The secular equation (3.14) then takes the form
3.18For an isotropic elastic half-space, we have
and then the secular equation reduces to the familiar form
3.19

## 4. The second-order problem

We now assume that the *x*_{3}=0 plane is a plane of material symmetry, which is the case for monoclinic materials. In this case, we may set . Also, there is no loss of generality in choosing and seeking a solution for in the form
4.1where *w*_{1} and *w*_{2}, both functions of *x*_{2} only, satisfy the following equations derived from (2.16) and (2.17):
4.2
4.3
4.4
and
4.5The four vectors **g**^{(1)},**g**^{(2)},**g**^{(3)} and **g**^{(4)} in the above equations are defined by their components
4.6To solve these two problems explicitly, we rewrite (3.7) in the form
4.7where the constants *b*_{αβ} can in principle be obtained by comparing (4.7) with (3.2) and (3.7). For instance, when the half-space is isotropic, these constants are given, to within a multiplicative constant, by
To solve the problem for *w*_{1}, we first note that a particular integral of (4.2) is given by
4.8where *s*_{1} and *s*_{2} are constants to be determined.

On substituting this into (4.2), making use of (4.7), and then comparing the coefficients of e^{ipβx2}, we obtain
4.9When the material is isotropic, the denominator of (4.9) becomes zero when *β*=1 owing to the fact that *p*_{1}=*p*_{3}. In this case, (4.8) is modified to
4.10however, for brevity, the corresponding expression for *s*_{1} is not written out here. Numerically, the case of isotropy can be dealt with by first considering a nearly isotropic material with e.g. moduli given by
4.11and then taking the limit *γ*→0.

A general solution of (4.2) is then given by
4.12where *p*_{3} is the root of
with positive imaginary part, and *s*_{3} is a disposable constant. We observe that such a *p*_{3} exists only if , that is if , where *v*_{t} is the anti-plane shear body wave speed. In other words, the type of topography-guided surface-wave solution under consideration only exists if the original surface wave on the associated flat half-space is subsonic. For the composites defined by (2.5) and (2.6) with , figure 2 shows the variation of and *v*_{t} with respect to the angle *θ*. It is seen that the subsonic condition is satisfied only for 0^{°}≤*θ*<24.37^{°} or 54.67^{°}<*θ*≤90^{°}.

In figure 3, we have shown the variation of and *v*_{t} with respect to the angle *θ* between the [100] direction and that of wave propagation, assuming that the surface waves propagate on the (001) plane of the silicon material described by (2.7). In this case, the surface wave is subsonic for all the values of *θ* considered [23].

On substituting (4.12) into the boundary condition (4.3), we obtain
4.13Similarly, the problem for *w*_{2} can be solved to yield
4.14where
4.15and
4.16Thus, under the assumption that *x*_{3}=0 is a plane of material symmetry, the second-order problem can be solved without having to impose any solvability conditions.

## 5. The third-order problem

If we write **u**^{(2)} as
5.1then the unknown vector function * y*(

*x*

_{2},

*x*

_{3}) satisfies the following equations derived from (2.18) and (2.19): 5.2and 5.3where a prime denotes partial differentiation with respect to

*x*

_{2}, 5.4 5.5and

*S*is the matrix with components

*c*

_{α3β3}.

It can easily be shown, by integrating by parts followed by the use of (3.3), that Thus, the solvability condition for the inhomogeneous problem (5.2) and (5.3) is given by 5.6We note that a similar solvability condition can be written down for the second-order problem, but it can be shown that it is automatically satisfied.

On substituting (5.4) and (5.5) into (5.6), we obtain
5.7where
5.8
5.9
5.10
5.11
and
5.12The coefficient *c*_{0} is obviously real and positive. In appendix A, we show that *c*_{2} and *c*_{4} are real, but *c*_{3} is purely imaginary and is related to *c*_{1} by *c*_{3}=2i Im(*c*_{1}). Thus, on writing *f* in the polar form
5.13it is then straightforward to deduce that and that the variation of *r*(*x*_{3}) is governed by the equation
5.14where
and and denote the real and imaginary parts of *c*_{1}, respectively.

Thus, the problem of finding the speed correction *X*_{1} is reduced to solving the eigenvalue problem of (5.14) subject to the decay conditions *r*(*x*_{3})→0 as . We recall, however, that the above eigenvalue problem has been derived under the assumption that the wavenumber *k* in (3.1) is unity. For the case when *k* is not unity, we may observe that under the substitutions *kx*_{1}→*x*_{1}, *kx*_{2}→*x*_{2} and *kx*_{3}→*x*_{3}, equations (2.14)–(2.19) remain the same, except that *h*′(*x*_{3}) and *h*′′(*x*_{3}) should be replaced by *h*′(*k*^{−1}*x*_{3}) and *k*^{−1}*h*′′(*k*^{−1}*x*_{3}), respectively. Equation (5.14) is then replaced by
5.15or equivalently, with another substitution *x*_{3}/*k*→*x*_{3} followed by *r*(*kx*_{3})→*r*(*x*_{3}),
5.16which shows that the speed correction is dependent on the wavenumber.

The amplitude equation (5.16) is recognized as a special case of the so-called time-independent Schrödinger equation
5.17associated with the potential well *V* (*x*_{3}) and energy level *E* (a constant, which should not be confused with the matrix *E* in (3.2)). Under the assumption that
Simon [25] and Klaus [26] have established the following general results concerning its negative eigenvalues.

(i) If there exists an eigenvalue when the potential well

*V*(*x*_{3}) is of sufficiently small amplitude, then the eigenvalue is unique and is given by 5.18where the star denotes convolution,*I*is defined by and*δ*is a small positive parameter characterizing the amplitude of*V*(*x*_{3}).(ii) A necessary and sufficient condition for the existence of the above-mentioned single eigenvalue is that the right-hand side of (5.18) is positive. Thus, under the assumption that the

*second term on the right-hand side of (5.18) is one order of magnitude smaller than the first term*, this condition is*I*≥0 (note that if*I*=0, the second term on the right-hand side of (5.18) is automatically positive; see Simon [25], p. 284).(iii) Since increasing the amplitude of a potential can only increase the number of eigenvalues [25], p. 284, the condition

*I*≥0 is also sufficient for the existence of at least one eigenvalue, even when the amplitude of*V*is not small.

For the special form of *V* in (5.16), we have
and so *I*≥0 if and only if *d*_{2}≥0. However, when *h* is of small amplitude, the two terms in the asymptotic expansion (5.18) are of the same order of magnitude, and they combine to give
5.19Thus, the existence condition is in fact given by
5.20This special case of *V* serves to demonstrate the fact that eigenvalues may still exist even if *I*<0 (whether *V* has small amplitude or not). We highlight this fact since some authors have previously claimed that *I*≥0 is also necessary for the existence of an eigenvalue. This claim is only valid if the condition in italics in item (ii) above is satisfied and if the potential is of sufficiently small amplitude.

For the case corresponding to figure 2, figures 4–6 show the variations of the three coefficients for values of the angle *θ* for which the two-dimensional surface wave is subsonic. The blow-up behaviour corresponds to the fact that *c*_{4} vanishes at *θ*=14.8^{°}*or*55.7^{°} approximately. It is seen that *d*_{2} is positive in both subsonic ranges (except at the two isolated values of *θ* where it blows up), and so it follows immediately from the above general results that (5.14) has at least one eigenvalue. This is confirmed in our numerical calculations later. In contrast, for the silicon material defined by equation (2.7), figure 7 shows that is negative for all values of *θ*, and so the existence of eigenvalues cannot be established using the general results above.

## 6. Isotropic materials

Before solving (5.16) subject to the decay conditions numerically, we first consider the special case when the material is isotropic. This case has previously been studied by Adams *et al.* [18] using a procedure that is particularly developed for isotropic materials. We now show that our formulae recover their results in this special case.

First, to facilitate comparison, we write *X*_{0} for *X*_{0}/*μ* throughout this section and introduce *κ* through . Equation (3.19) can be rewritten as
6.1The second-order problem now reduces to
6.2
6.3
6.4
and
6.5where
It is then easy to show that the solutions are given by
with
As a result, the expressions (5.8)–(5.12) can be evaluated explicitly to give *c*_{3}=0, *c*_{2}=0,
6.6where
and
By cross-multiplication followed by repeated use of (6.1) to eliminate powers of *X*_{0} higher than 3, we have verified using Mathematica that our *c*_{4}/*c*_{1} and *c*_{4}/*c*_{0} are equal to those of Adams *et al.* [18], *βA*/*B* and *A*/*C*, respectively. Additionally, as a check on our derivations, we have used the elastic moduli given by (4.11) to compute the coefficients for increasingly small *γ* and have obtained the same values as the above explicit expressions.

We observe that with *d*_{2}=0, the sufficient condition (5.20) is satisfied automatically, and so we may conclude that topography-guided surface waves always exist in an isotropic material. This settles the existence question left open in Adams *et al.* [18]. We may further conclude that when the geometric inhomogeneity is of sufficiently small amplitude, there can only exist a single guided surface wave. Finally, it can easily be verified numerically that *c*_{0}/*c*_{4} (=*d*_{0}) is positive for all . We thus conclude that all topography-guided surface waves on an isotropic elastic half-space travel at a slower speed than the classical Rayleigh wave (see (7.2) and also Bonnet-Ben *et al.* [12], Theorem 1).

## 7. Numerical results

We shall now explain our numerical procedure by focusing on the case *k*=1. Thus, we return to the amplitude equation (5.14) and solve it subject to the decay conditions . In the limit , equation (5.14) can be approximated by
7.1It is clear that *r*(*x*_{3}) will have the required decay behaviour as only if
7.2and when this is satisfied, we have
7.3It then follows that if *d*_{0}<0, then the topography-guided surface waves travel at a lower speed than their two-dimensional counterpart, whereas if *d*_{0}>0, then the topography-guided surface waves are faster than their two-dimensional counterpart. When *h*(*x*_{3}) is an even function of *x*_{3}, *r*(−*x*_{3}) is a solution of (5.7) whenever *r*(*x*_{3}) is a solution. Thus, the eigensolutions of (5.14) are either even or odd. For the even (symmetric) modes, we may impose, without loss of generality, the conditions
7.4and the decay behaviour through
7.5where *L* is a sufficiently large positive constant and the first equation in (7.5) defines the error function *e*(*X*_{1}). For each fixed *X*_{1}, the *e*(*X*_{1}) can be evaluated after integrating (5.14) subjected to the initial conditions (7.4). We first plot *e*(*X*_{1}) against *X*_{1} to show the approximate locations of any zeros and then use the Newton–Raphson method to find the exact values of the zeros.

For the odd (anti-symmetric) modes, the initial conditions (7.4) are replaced by
7.6and we integrate (5.14) subject to the initial conditions (7.6) and iterate on *X*_{1} in order to satisfy the decay condition (7.5).

When *h*(*x*_{3}) is not an even function of *x*_{3}, the eigenmodes do not have any symmetry properties. We may, for instance, integrate (5.14) subject to (7.5) and
7.7respectively, and iterate on *X*_{1} so that the two solutions have the same gradient at a matching point, say *x*_{3}=0.

The above numerical scheme is first tested on the eigenvalue problem
7.8which is known to have *n* eigenvalues given by *λ*=*i*^{2} (*i*=1,2,…,*n*), with odd *i* corresponding to odd modes and even *i* to even modes; see Drazin & Johnson [27]. Our scheme is able to reproduce these exact results correctly. The numerical scheme is then applied to solve (5.14) for a variety of *h*(*x*_{3}) and different values of *d*_{1} and *d*_{2}. Our numerical results confirm the validity of the asymptotic formula (5.19) and the existence condition (5.20). In particular, we invariably find that for *h* sufficiently small, the single eigenvalue tends to zero when approaches zero from above, and there is no eigenvalue when .

We now present illustrative calculations for the composite material defined by (2.6) and for the case when the topography is described by the ‘Gaussian bump’ . When *θ*=2.187^{°}, we have *d*_{0}=0.36 and the corresponding *X*_{1} must necessarily be negative. It is found that there is only a single solution that is symmetric with *X*_{1}=−0.146 (the unit of *X*_{1} is 10^{9} *N* m^{−2}), and so the topography has the effect of reducing the two-dimensional surface-wave speed. The corresponding profile of *r*(*x*_{3}) is shown in figure 8.

When *θ*=16.767^{°}, we have *d*_{0}=−0.77, and the corresponding *X*_{1} must necessarily be positive. It is found that there exist one odd solution and one even solution, corresponding to *X*_{1}=0.186 and 0.469, respectively. The topography has the effect of producing two variations of the original two-dimensional surface wave. Both waves have a higher speed than their two-dimensional counterpart, and their lateral variations are shown in figure 9.

We have also carried out calculations for the silicon material defined by (2.7) for which for all values of *θ*. We have tried a large number of *h* profiles, but have found no eigenvalues for (5.14).

## 8. Concluding remarks

We have extended the analysis of Adams *et al.* [18] to a generally anisotropic elastic material with *x*_{3}=0 a plane of material symmetry. It is shown that determination of the wave speed correction is also reduced to the solution of a simple eigenvalue problem. However, the current anisotropic problem has a much richer variety of behaviour and raises new questions. For instance, the existence of a topography-guided surface wave can be eliminated if the original two-dimensional surface wave is supersonic with respect to anti-plane shear motions, and surface inhomogeneity may either increase or decrease the speed of the original two-dimensional surface wave. There also exists the possibility that the coefficient of the second-order derivative term in the reduced eigenvalue problem vanishes, in which case, a topography-guided surface wave, if it exists, may be expected to behave very differently.

The methodology that we employed is based on the assumption that a slowly varying perturbation to the otherwise flat surface will give rise to a small-amplitude perturbation to the surface-wave speed. This may, however, not be the case. For instance, Farnell [23] showed that as the direction of propagation is slightly turned in the *x*_{1}*x*_{3}-plane (in terms of the notation in the present paper), a two-dimensional surface wave will become a three-dimensional surface wave, and the wave speed may experience a finite change, no matter how small the turning is. Thus, when the associated two-dimensional surface wave is supersonic with respect to anti-plane shear motions, we cannot exclude the existence of a topography-guided surface wave whose wave speed differs significantly from the two-dimensional surface-wave speed. Finally, when the *x*_{3}=0 plane is not a plane of material symmetry, the present methodology does not apply either. We leave these questions to future studies.

## Acknowledgements

The contribution of W.F.W. was supported by a PhD studentship jointly funded by China Scholarship Council and Keele University. We thank all three referees for their constructive comments and suggestions.

## Appendix A. Proof that *c*_{2} and *c*_{4} are real and *c*_{3}=2i Im (*c*_{1})

We first define a differential operator by
A1see (4.2). Then for any twice differentiable function *v*(*x*_{2}), we have, by integrating by parts,
A2where an overline signifies complex conjugation. With *v* replaced by *w*_{1}, the above identity shows that the expression
is real. This implies, through the further use of (4.2) and (4.3), that the expression
is real. The reality of *c*_{4} then follows from this result and the fact that *S* is a symmetric matrix so that the extra term is also real.

In a similar manner, upon replacing *v* and *w*_{1} in (A2) both by *w*_{2} and making use of (4.4) and (4.5), we find that *c*_{2} is also real.

If, on the other hand, the *v* in (A2) is replaced by *w*_{2} and use is made of (4.2)–(4.5), we obtain
A3The relation then follows from this result after straightforward manipulations.

- Received June 21, 2012.
- Accepted September 25, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.