A Factorizing {\cal K}(z)

Dynamic interfacial trapping of flexural waves in structured plates


The paper presents new results on localisation and transmission of flexural waves in a structured plate containing a semi-infinite two-dimensional array of rigid pins. In particular, surface waves are identified and studied at the interface boundary between the homogeneous part of the flexural plate and the part occupied by rigid pins. A formal connection has been made with the dispersion properties of flexural Bloch waves in an infinite doubly periodic array of rigid pins. Special attention is given to regimes corresponding to standing waves of different types as well as Dirac-like points, that may occur on the dispersion surfaces. A single half-grating problem, hitherto unreported in the literature, is also shown to bring interesting solutions.

Department of Mathematical Sciences, University of Liverpool, Peach Street, Liverpool L69 7ZL, United Kingdom

Department of Mathematics, Imperial College London, London SW7 2AZ, UK

School of Engineering, Liverpool John Moores University, Liverpool L3 3AF, UK

1 Introduction

The advent of designer materials such as metamaterials, photonic crystals and micro-structured media that are able to generate effects unobtainable by natural media, such as Pendry’s flat lens [24], is driving a revolution in materials science. Many of these ideas originate in electromagnetism and optics, but are now percolating into other wave systems such as those of elasticity, acoustics or the idealised Kirchhoff-Love plate equations for flexural waves, with this analogue of photonic crystals being labelled as platonics [18]. Many of the effects from photonics also appear in the flexural wave context albeit with some changes due to the biharmonic nature of the Kirchhoff-Love equation: ultra-refraction and negative refraction [6, 7], Dirac-like cones, Dirac-cone cloaking and related effects [3, 27, 19] amongst others.

The Kirchhoff-Love equations are good approximations within their realm of applicability (for instance with , as plate thickness and typical wavelength respectively then is required) and capture much of the essence of the wave physics, hence their emergent popularity. The system is relatively simple so analytic results for infinite periodic structured plates pinned, say, at regular points readily emerge [16], with much of this earlier work reviewed in [20]; more recently multipole methods [21], extending pins to cylinders, or high-frequency homogenisation approaches [2] to get effective continuum equations that encapsulate the microstructure, have emerged. For finite pinned regions of a plate, a Green’s function approach [5] leads to rapid numerical solutions, or for an infinite grating, one may employ an elegant methodology for exploring Rayleigh-Bloch modes. This includes extensions to stacks of gratings and the trapping and filtering of waves [22, 25, 11] which further exemplify this approach. Problems such as these start to pose questions about semi-infinite gratings, or edge states in semi-infinite lattices, and our aim is to generate the relevant exact solutions.

For semi-infinite cracks, and other situations where the change in boundary condition along a line is of interest, the classical Wiener-Hopf technique [31] for solving integral equations is highly developed and often used for mixed boundary value problems in continuum mechanics; indeed the seminal Wiener-Hopf treatise by [23] covers this aspect almost exclusively. However, lesser known and presented only as an exercise (4.10, p.173-4 [8]) in [23], is its application to continuum discrete problems such as gratings. As such, it has seen notable application to the fracture of discrete lattice systems as reviewed by [26], to antenna design [30] and various semi-infinite grating and lattice scattering problems in acoustics [15, 29] where the early work of [12] (hereafter referred to as HK), motivated by [14], provides a wealth of useful information.

The corresponding exact solutions for the platonic system are not available and, given the current interest in platonics and the versatility of these in, say, asymptotic schemes and the key insight given into the physics and results, we aim to provide these here. Of particular interest are regimes for which one observes the interplay between grating modes and Floquet-Bloch waves in a full doubly periodic infinite structure. In particular, these include frequencies and wave vector components corresponding to stationary points on the dispersion surfaces as well as neighbourhoods of the Dirac-like points. Additionally, we address a fundamental question as to whether a simple, yet surprisingly physically rich, structure such as a half-plane of rigid pins in a plate supports flexural interfacial modes; by which we mean waves that propagate along the interface of the platonic crystal and the homogeneous part of the biharmonic plate.

Figure 1: (a) A single semi-infinite grating of rigid pins in a biharmonic plate with spacing . The pins are represented by the large black discs. The forces/intensities are represented by the coefficients , and displacements, used in the Wiener-Hopf approach, by . (b) The semi-infinite lattice of pins, consisting of an array of infinite gratings with period in the vertical direction and horizontal spacing .

Since the underlying mathematical structure of a semi-infinite linear grating and that of a semi-infinite lattice are very similar we choose to treat both examples within this paper. For clarity, we shall refer to a grating when we have a single half-line of scatterers (see figure 1(a)) and to a lattice when there is a lattice of scatterers in a half-plane (see figure 1(b)); throughout this paper we shall assume that we have point clamped scatterers although the analysis is readily extended to more complicated conditions holding at a point. In section 2, we formulate the problem drawing upon the discrete Wiener-Hopf technique used by HK for the semi-infinite diffraction grating governed by a Helmholtz operator; much is directly applicable provided we replace the point source Green’s function for the two-dimensional Helmholtz operator with that for the biharmonic operator. In some regards, difficulties associated with the Helmholtz Green’s function’s logarithmic singularity are bypassed for the biharmonic case, where there is no singularity. We turn our attention to the linear grating first, in section 2.1, and follow that with an analysis of the lattice in section 2.2, as noted above this is a parametrically and physically rich problem displaying a surprisingly wide range of behaviours and we outline methodologies for both interpreting and designing specific behaviours; we consider these in sections 3.13.2 for the semi-infinite grating and lattice respectively. Concluding remarks are drawn together in section 4.

2 Formulation

The methodology primarily adopted here, at least for the exact solution, is that of discrete Wiener-Hopf drawing upon the treatment of the Helmholtz semi-infinite grating by HK, which we briefly review here. To avoid issues with singularities in their problem, HK consider cylinders of finite radius, small compared to wavelength, and widely spaced. The diffracted field consists of a cylindrical wave plus a set of plane waves that possess amplitudes and directions of propagation identical to those for an infinite grating, but for the semi-infinite case, they do not exist everywhere. Instead the end-effects lead to shadow boundaries defined by lines drawn from the end of the grating to infinity along propagation directions of various plane waves. The cylindrical wave has sharp zeros followed by a maximum for certain directions; these zeros occur in the directions that would be shadow boundaries if the incident wave travels directly into the end of the grating - a class of resonant cases characterised by a diffracted wave that travels parallel to the grating. HK note that for the wave that travels parallel and out of the grating, the cylindrical wave and all the other plane waves vanish.

Assuming isotropic scatterers HK, see also [15], pose the solution in the form


for some coefficients to be determined and being the usual Hankel function of the first kind with wavenumber , and representing the distance between the observation point and the th scattering element on the -axis; for an infinite array these just differ by a phase factor and the analysis is simpler. The assumption that the grating’s elements are small, and therefore scatter isotropically, are not formally required for the pinned biharmonic plate as this is exactly true.

In HK, the spacing between the scatterers is assumed to be large compared with the wavelength, and it must not be an integral or half-integral multiple of the wavelength. These assumptions are required by the mathematical analysis employed. The reason for the limitation on the integral or half-integral multiples of wavelength is to ensure that branch cuts for and are distinct (see appendix B for a discussion of branch cuts for the biharmonic plate) and is the spacing of the pins.

To determine the scattered wave coefficients in (1), the boundary conditions generate an infinite system of linear equations solvable using discrete Wiener-Hopf. In the neighbourhood of the th scatterer, the field consists of the incident plane wave, the sum of waves scattered to the th element by the other elements and the wave scattered by the th element itself:


Here is the radius of the individual elements. The term is treated separately because of the logarithmic singularity arising for the Hankel function at the origin. Although a small radius is assumed by HK, the singularity poses problems as and much of HK involves dealing with this; for the biharmonic case we discuss in this article, the Green’s function is finite at the source point rather than diverging logarithmically.

2.1 Semi-infinite grating of rigid pins in a biharmonic plate

For the Kirchhoff-Love plate equation


with an array of pins, the solution, cf. (1), is


with the to be determined, and where the free-space Green’s function


is a function of position , with for spacing and . As in (1) the Green’s function contains but this is now augmented by which is the modified Bessel function [1]; this Green’s function is bounded at .

We take rigid pins located at points for and to set up the system of equations we introduce displacements at for all ; imposing zero displacements at the pins, for but are unknown for . For the intensities in (4) since there are no sources on the left-hand side, . The analogue of (2) is


We now employ the -transform, multiply (6) by for complex, and sum over all :


To transform (7) into a single functional equation of the Wiener-Hopf type, it is convenient to define and a kernel function as


Notably, for , the term is constant, , which replaces the awkward term of (2). Using these functions, (7) becomes


and the forcing function, in this case, is


but could take different forms if the forcing were altered.

Equation (10) is the starting point for the Wiener-Hopf technique with key steps being the product factorization of the kernel function into two factors which are analytic in given but different regions, and factorization of the forcing function. Before proceeding we turn to the semi-infinite lattice as the formulation is almost identical.

2.2 Semi-infinite lattice of rigid pins in a biharmonic plate

We now replace each pin in the semi-infinite grating with an infinite grating in the vertical direction, as illustrated in figure 1(b). The properties of the infinite grating are well known e.g. [10] and the quasi-periodic grating Green’s function


plays an important role. Here is the vertical period (in section 2.1 we used to denote the period for the single grating) and is the corresponding Bloch parameter. The gratings are separated by integer multiples of spacing and the semi-infinite lattice is created from columns of infinite gratings (or rows of semi-infinite gratings). The analogue of (6) emerges as


We repeat the procedure for the semi-infinite grating with the only difference being that the kernel function is now connected with the doubly quasi-periodic Green’s function:


when with .

Figure 2: The unit circle, , and contours and which define the boundary of the ring of analyticity. The integrals on and are evaluated counter-clockwise and clockwise respectively.

2.3 Wiener-Hopf

The crux of the Wiener-Hopf methodology is to unravel (10) using the analyticity properties of the unknown functions and the known functions , . To proceed we define domains and such that


where is a circle of radius and is a circle of radius for (see figure 2); the intersection describes an annulus of analyticity in the neighbourhood of the unit circle. From the definitions of (8) and (10), and assuming appropriate decay at infinity, we identify them as and functions analytic in and respectively. We then attempt to separate (10) into sides entirely analytic in ; the only way both can be equal, given an overlapping region of analyticity, is for them both to equal the same entire function - thus identifying and . The outcome is


which requires the product factorization and sum factorization where ; the product factorization is highly technical and is relegated to appendix A, the sum factorization is straightforward by inspection. After some algebra we obtain


except for different definitions of , and , this mirrors HK. The common entire function is identified, after using Liouville’s Theorem to extend to the whole complex plane, as zero. Thence


and the -transform for the displacement coefficients also follows.

Technical details of the discrete Wiener-Hopf method used here are given in the appendices. The kernel function is the sum of an infinite series of Hankel and modified Bessel functions, and its factorization is discussed in appendix A, together with details about the regularisation techniques used for evaluating various integrals. Appendix B provides explanations and formulae required to accelerate the extremely slow convergence of the kernel’s series, owing to the highly oscillatory nature of the Hankel function terms and the presence of branch cuts. We also suggest some alternative approaches to accelerate numerical evaluation of the series besides the regularisation method primarily implemented here.

3 Results

We use two approaches to present results and illustrative examples for the displacement fields associated with the two types of array; in section 3.1 we consider the single semi-infinite pinned platonic grating characterised by spacing , and in 3.2 we analyse the two-dimensional lattice defined by and . One of the approaches is the discrete Wiener-Hopf method described above and we compare the results with those for a truncated semi-infinite array analysed with a method attributable to [9], which we outline in section 3.1 for the single grating.

3.1 The semi-infinite grating of rigid pins

For a plane wave incident at an angle as in figure 1(a), we determine the coefficients for a truncated grating by solving the algebraic system of linear equations:


with being the spacing of the grating’s pinned points and the incident wave as defined in equation (4). This is the standard Foldy scattering equation ([9]), and used amongst others by [15] for the Helmholtz problem, and [5] for the biharmonic plate. It is solved for pinned points, and the displacements are plotted using


We compare results with those obtained from the discrete Wiener-Hopf technique for the semi-infinite grating, which gives us the exact solution. The expressions for  (A3) are substituted into equation (18) to determine , bearing in mind the additional singularity arising from the term in the denominator. Referring to the expansion (8) we take the inverse of the -transform to determine the coefficients . Multiplying (8) through by and integrating with respect to , we obtain


The final result for the coefficients is the integral


for which the singularities and the branch cuts of appendix B have to be taken into account. The results compare well with a truncated version (at least 2000 pins) treated with the Foldy scattering approach. In the examples that follow we use both real parts and moduli of the displacement field defined by these complex coefficients, since they give us insight into the propagation tendencies of the scattered waves.

For the far-field behaviour of the coefficients we deduce a relation of the form


The case denotes the propagating Bloch wave and for we observe localisation linked with the exponential decay of the coefficients. The ratio is determined by


where the kernel in its present form is evaluated numerically.

Resonant cases

We consider some of the frequency regimes mentioned by [5] for a finite array of pinned points and an infinite grating, as well as those referred to by HK for a semi-infinite grating, including the case they define as a resonance. This is where the incident wave travels directly into the end of the grating ( for instance) and then the diffracted wave travels parallel to the grating, either into (inward resonance) or away from the grating (outward resonance).

Spectral orders of diffraction are defined according to the equation (see for instance [4])


with only a finite number of the being real and representing propagating waves. The remaining orders are complex and represent evanescent waves. According to HK, the resonant cases arise when one of the spectral directions is zero (inward resonance) or , the case of outward resonance. Thus, resonances coincide with additional diffraction orders becoming propagating. For outward resonance


thus for , and outward resonances occur for frequencies corresponding to being a multiple of .

Figure 3: Plane wave incident on a semi-infinite grating with spacing at . (a-c) Real part of the scattered displacement , for . (d) Comparison of the moduli of the coefficients using Wiener-Hopf (dots for ) and Foldy methods (shown for all).

We observe some evidence of this effect for the case in figure 3, where we plot the scattered displacement field for values of around for angle of incidence and period . We use Foldy truncated to pins and the results from Wiener-Hopf and Foldy are visually indistinguishable for the first several gratings, as shown in figure 3(d) which compares the moduli of the coefficients in both cases for (i.e. the black solid line and dots for Foldy and Wiener-Hopf respectively).

In figure 3(a-c) we show the real part of scattered displacement , for and , that is, we pass through the resonance at . For shown in figure 3(a), the high reflected energy is illustrative of the outward resonance described by HK for the resonant value .

We contrast this with in figure 3(c) for which an additional diffraction order has become propagating. For diffraction orders to pass off, . Thus for , the order becomes propagating along with for the resonant value , and the presence of two distinct propagating orders is clearly illustrated for . The two examples either side of show strong evidence of the circular wave that is typical for the end-effects of a semi-infinite scatterer.

The resonant frequency case of is shown in figure 3(b) and there are relatively lower amplitudes of displacement (with a maximum of ) with most of the scattering occurring along the grating, contrasting with the effects illustrated in figures 3(a) and (c). The localisation and reduction of the maximum amplitudes to is consistent with HK’s observation that for outward resonance, the cylindrical wave and all but one of the plane waves vanish. Indeed, the real part of the total displacement field (not shown) closely resembles that of the incident wave, except for the region surrounding the grating itself. Comparing the coefficients , as we do in figure 3(d), further emphasises the difference between this resonant case, for which rapidly saturates to a low value versus those for that take much longer to saturate and then do so to a larger value.

Shadow boundaries

HK refer to non-resonant cases where the diffracted field for a Helmholtz-governed semi-infinite grating consists of a set of plane waves and a cylindrical wave. The plane waves are consistent with those arising for the infinite grating, but do not exist everywhere. We obtain similar results for the biharmonic case, for which the propagating waves are due to the Hankel functions arising from the Helmholtz part of the Green’s function in equation (4). The coefficients for the scattered field are defined by equations (18) and (22) leading to the expression:


Furthermore, the representation (4) for the displacement, leads to the following approximate expression for large values of :

HK have proved in their paper that the shadow boundaries correspond to singularities of the -transform (18), when is a point on the unit circle defined in the form . We note that the regularised kernel does not have roots and poles on the unit circle, and hence the shadow boundaries are defined by the equation (25) and coincide with those discussed by HK in [12].

Figure 4: The scattering of a plane wave by a grating with , for . (a) Real part of displacement field . The shadow boundaries for and are marked using dashed black lines. (b) Comparison of Wiener-Hopf (dots) and Foldy (solid black) coefficients .

Namely, denoting , there are poles whenever . The resulting values of are labelled by and are the directions of the spectral orders of an infinite grating (25), of which a finite number are real and propagating. In computing the residues at the poles, it is merely sufficient to note that the factorization exists rather than evaluating the factors explicitly.

The lines act as shadow boundaries and this is illustrated in figure 4 where and . For , the diffraction order becomes propagating for , so for , both and are propagating. We observe the two shadow boundaries for and , and the field consisting of two propagating diffracted waves, as well as the reflected wave. Figure 4(b) complements figure 4(a) by showing a comparison of the moduli of the coefficients with the Foldy method represented by the solid black line, the Wiener-Hopf by the dots; for oblique incidence, the two methods match well and for the Wiener-Hopf we use the regularisation parameter and 1200 intervals for the trapezoidal integration (see appendices AB). Increasing the number of intervals leads to closer convergence to the Foldy coefficients.

Reflection and transmission

Frequency regimes associated with total reflection and total transmission by an infinite grating, as identified by [5, 22], are interesting cases. For an infinite pinned grating (period , ) the normalised reflected and transmitted energies are plotted against in figure 5(a).


Figure 5: Plane wave incident on a semi-infinite grating with spacing at an angle . (a) The normalised reflected, , (solid blue curve) and transmitted, , (dashed black) energy versus for scattering by an infinite pinned grating. The real part of the displacement field for the semi-infinite grating is shown for (b) , (c) and (d) .

Three important values of = 3.2; 3.55; 3.68, which are total reflection, equipartition of energy and total transmission respectively, are labelled b-d and we show the real part of the total displacement field for the corresponding semi-infinite gratings in figures 5(b)-(d).

Figure 5(a) shows that the maximum reflected energy (solid blue line) for the zeroth propagating order arises for . The Wood anomaly at signifies the passing off of the order , explaining the multiple orders illustrated in figure 4 for . The real part of the total displacement field for the semi-infinite grating for the same parameter values and shows strong reflection to the right of the vertex, consistent with that observed for the full grating for the single propagating plane wave. This contrasts sharply with the results observed for the frequency associated with in figure 4. Note that there is only one shadow boundary in figure 5(b) (for the shadow region behind the grating) whereas figure 4(a) shows evidence of two distinct shadow lines.

In figure 5(c) we highlight which supports a mixture of reflected and transmitted energy for both the infinite and semi-infinite gratings; the dashed line in figure 5(a) represents the transmitted energy. The real part of the displacement field in figure 5(c) is consistent with the information provided by the energy diagrams for the corresponding infinite grating; both reflection and transmission are visible, with the shadow region above the grating now admitting transmittance. This transmission effect dilutes the reflection of figure 5(b), and total reflectance is converted to total transmission by adjusting the parameter appropriately. This is illustrated in figure 5(a) where leads to full transmission for the infinite grating, and produces a similar effect for the semi-infinite grating whose total displacement field is shown in figure 5(d). This transmission resonance is an example of the outward resonant effect for oblique incidence. The scattered field (not shown) is reminiscent of figure 3(b); evidence of the outward-travelling wave close to the vertex is visible in figure 5(c) for .

3.2 Semi-infinite lattice of rigid pins

For an incident plane wave characterised by as in figure 1(a), we determine the coefficients for a truncated half-plane by solving the algebraic system of linear equations (19) but replacing each pin with a grating in the vertical direction. We compare results with those obtained using discrete Wiener-Hopf.

Zeros of the kernel function

Recall from equations (12), (14) that for lying on the unit circle, the function is precisely the doubly quasi-periodic Green’s function. Therefore its zeros represent points on the dispersion surfaces for Bloch waves in the infinite doubly periodic structure. This direct connection between the kernel and the doubly periodic medium enables us to analyse wave phenomena using dispersion surfaces, band diagrams and slowness contours.

We determine the scattering coefficients and displacements using the discrete Wiener-Hopf method. A good approximation to these results is obtained by applying Foldy’s method to a large enough array. For systems of at least gratings we demonstrate several interesting wave phenomena. [19] highlight the potential for varying the aspect ratio of rectangular lattices of pins, linking Dirac-like points with parabolic profiles in their neighbourhood. The characteristic of this parabolic profile determines the direction of propagation of localised waves. We consider first a semi-infinite rectangular lattice with aspect ratio . The band diagram for the doubly periodic rectangular array is shown in figure 7 of [19]. Here we represent the dispersion surfaces with isofrequency contour diagrams.

For our initial investigations we focus on the first three dispersion surfaces.

Figure 6: Isofrequency contours (a selection of values are labelled) for the third dispersion surface of the rectangular lattice with , . The Brillouin zone is marked by the solid blue rectangle, with .

In figure 6 we illustrate the third surface for the rectangular lattice with and with the Brillouin zone labelled by , where . We observe a parabolic cylinder profile along with an inflexion at where the contours change direction for and a Dirac-like point for at . We also note that the concentration of points in figure 6 tending towards the straight lines

and periodic shifts by , are projections (onto the -plane) of intersections of two light cones, which are given by


The remaining surfaces are presented together with displacement fields in the examples that follow. Both for the first surface and for the second surface are mentioned by [19]. This immediately gives us some frequency regimes to investigate for finding wave phenomena including neutrality and interfacial waves. We also demonstrate blocking and resonant behaviour.

Wave-vector diagrams

[32] demonstrated an elegant strategy using wave-vector diagrams to investigate wave phenomena in planar waveguides; the technique was also outlined by [13] for photonic crystals in their Chapter 10. [32] was able to illustrate negative refraction, focusing and interference effects, having predicted them with careful analysis of the wave vector diagrams. We adopt an analogous approach for this platonic crystal system.

The underlying principle is a corollary of Bloch’s theorem; in a linear system with discrete translational symmetry the Bloch wave vector is conserved as waves propagate, up to the addition of reciprocal lattice vectors. Since there is only translational symmetry along directions parallel to the interface, only the wave vector parallel to the interface, , is conserved. In our case, this is precisely the direction . Thus for any incident plane wave defined by and , any reflected or refracted (transmitted) wave must also possess the same frequency and wave vector for any integer and some . Therefore, wave-vector diagrams may be used to analyse the reflection and refraction of waves within the pinned system.

These isofrequency diagrams (also called slowness contour plots) consist of dispersion curves for constant (characterising the platonic crystal) and the contour for the ambient medium (the homogeneous biharmonic plate) on the same , diagram (see figure 6). The incident wave vector (whose group velocity direction is characterised by ) is appended to the ambient medium’s contour (the circle in figure 6), and a dashed line perpendicular to the interface of the platonic crystal and the homogeneous part of the biharmonic plate ( or direction here) is drawn through this point. The dashed line and have been added to figure 6. The places where this dashed line intersects the platonic crystal contours determine the refracted waves and their directions, such that their group velocity is perpendicular to the contours, and points in the direction of increasing .

Additional information relating to the coefficients is obtained from the Wiener-Hopf method. Exponential decay of the form (23) with would indicate the possibility of interfacial waveguide modes, with verification supported by the prediction of the refracted waves’ directions using wave-vector diagrams. We illustrate the technique with an introductory example demonstrating a localised mode along the edge of the crystal for the first dispersion surface, as well as reflection and transmission which is predicted from the isofrequency diagram.

Reflection, transmission and interfacial waves for the first dispersion surface

In figure 7(a) we show a collection of isofrequency contours for constant for the first dispersion surface for the rectangular lattice with aspect ratio .

Figure 7: (a) Isofrequency contours constant for the first dispersion surface for the rectangular lattice with , . (b) Isofrequency contours , 3.107, 3.11.

This range of frequencies is notable for its flat bands and the parabolic profile mentioned by [19] at in the vicinity of a Dirac-like point. For normal incidence , indicated by the direction of the arrow in figure 7(a), the dashed line intersecting the ambient medium’s contour and normal to the interface () corresponds to , as shown in figure 7(a). We seek the intersections of this line with the isofrequency contour for a specific choice of . Since the direction of the group velocity of the refracted waves is perpendicular to these constant contours, for an interfacial mode we would require the contour to be tangent to the line .

Interfacial waves

It is clear from figure 7(a) that the contours tend towards tangency as increases to but then their behaviour in the vicinity of is less predictable. This is connected both with the resonant case for the single semi-infinite grating, see section 3.1.1, and the parabolic profiles associated with Dirac-like points to which [19] alluded for . This observation is supported by the discontinuous nature of the contours for when they intersect the line in figure 7(a), and in particular the two distinct contour curves labelled by on the right. The contours change direction at for ; this corresponds to a point of inflexion on the band diagram.

In figure 7(b) we show isofrequency contours for (blue), 3.107 (orange) and 3.11 (red) to emphasise the transition of the contours through this narrow frequency window. The curve for is consistent with preceding values of but it does not touch the line, whereas (orange curve) is extremely close to touching, with the origin of a point of discontinuity clearly observable. This point of inflexion for some is the limit as the contours tend to and occurs at . For this value of , the refracted wave would travel along the interface, but there would be no preferential direction of propagation since the upper contour would indicate a “downward” direction (increasing ) whereas the lower contour would indicate the opposite direction. This suggests the presence of a standing wave, localised within the first few gratings of the half-plane of pins.

Figure 8: A plane wave is incident at on a lattice of 2000 gratings with , . (a-c) Real part of total displacement field for , , . (d) Comparison of moduli of coefficients .

We see evidence of this when we plot the real part of the displacement field for in figure 8(b), with the predicted preferential direction of propagation indicated in figure 7(b). The moduli of the coefficients are plotted versus -position of grating in figure 8(d), for the three values of featured in figure 7(b). Both and show exponential decay of the coefficients, suggesting localisation of waves, whereas the steady oscillation of the coefficients for predicts transmission of the waves through the system. Figure 8(b) clearly illustrates the localisation of waves within the first five gratings, and a clearer example of an interfacial wave is illustrated in part (a) for , although comparison with lower values of shows similar results, but with more striking reflection action.


Figure 7(b) predicts that for , the preferred direction of propagation of the refracted waves is normal to that for , in the form of transmission through the system. Note that the opposite direction ( rather than 0) is ruled out since any intersections corresponding to the group velocity directed towards the interface from the crystal would violate the boundary conditions. The coefficients’ behaviour in figure 8(d) also supports the hypothesis of propagating waves. This is demonstrated in figure 8(c) where we observe transmission as the wave propagates through the system without any change in direction, with the period of the wave’s envelope function consistent with the period for the coefficients in part (d).

Coupling with finite grating stacks

The spikes in figure 8(d) for and seem to indicate resonant interaction with the first few gratings of the semi-infinite array, as does the localisation evident in figures 8(a, b). There is a connection with the finite grating stacks analysed by [10]; resonances associated with the Bloch modes for the finite systems are linked to the neighbourhood of the Dirac-like point illustrated in figures 7(a, b).

For the pinned waveguide consisting of an odd number of gratings of period , a quasi-periodic Green’s function (12) is used to derive the dispersion equation for Bloch modes within the system. At each pin, the boundary conditions are . Therefore for a system of aligned gratings with spacing ,


where is the displacement, are coefficients to be determined for each Green’s function and is the number of gratings in the finite stack. This is equivalent to the matrix equation


where is the column vector of coefficients and the Green’s function matrix is complex and symmetric Toeplitz. The accompanying dispersion equation is:


the solutions of which characterise the system’s Bloch modes.

Figure 9: (a) Solutions of the eigenvalue problem for a system of five gratings with period and separation . (b) Close-up of frequency window for Bloch mode frequencies; local maximum and root marked.
Figure 10: A plane wave is incident at on a lattice of 2000 gratings with , . (a) Real part of coefficients for . (b) Real part of total displacement field for . (c) Real part of coefficients for . (d) Real part of total displacement field for .

Figure 8(b) features localised waves within the first five gratings so we solve the eigenvalue problem for a system of five gratings with period and separation in figure 9. The frequencies for the 5-grating system’s Bloch modes coincide with the frequency window for the Dirac-like point of the doubly periodic rectangular array. Similar results have been obtained for all -grating stacks with , with localised modes always confined to the range . This is consistent with figures 7(a, b) where the contours are well-behaved for and , but exhibit discontinuities in the neighbourhood of the Dirac-like point.

It has been shown that a finite system’s resonant Bloch modes may be coupled to appropriately chosen incident plane waves to generate transmission action [10]. We observe the same phenomenon for the half-plane in figure 8(c) as well as evidence of similar localisation for , in parts (a) and (b). Figure 9(b) adds further insight when we select two specific values of marked on the diagram; the local maximum for and the root which is close to the Dirac-like point value . In figure 10(a) the real parts of the coefficients for are of order , whereas for in figure 10(c), they are . We plot the accompanying total displacement fields in figure 10(b, d) where one should take note of the scaling bars. Both values in the vicinity of the Dirac-like point support propagation, but the intensities of the transmitted waves are linked to the turning points and roots for the finite grating stack structures, illustrated in figure 9.

This coupling can be switched off by altering the period and/or separation of the gratings or the incoming angle of incidence. It should be noted that these examples for the lowest frequency band involve the zeroth propagating order only. Higher frequencies introduce additional diffraction orders and their coupling facilitates more exotic wave phenomena including negative refraction, similar to those recorded by [32] for optical waveguides.

Neutrality in the vicinity of Dirac-like cones

We now consider neutrality effects that arise for parabolic profiles in the vicinity of Dirac-like cones. It is well known that close to Dirac points, waves may propagate as in free space, unaffected by any interaction with the microstructure within the crystal medium. These directions of neutral propagation around the Dirac point engender cloaking properties within the crystal, meaning that there is potential for “hiding” objects within the appropriate frequency regime. This property of neutrality for platonic crystals was mentioned by [18], and was explained in a simple way in terms of the singular directions of the Green’s functions for the biharmonic equation.

Figure 11: A plane wave is incident at on a lattice of 4000 gratings with , for . (a) Real parts of the coefficients . (b) Real part of the scattered field.

The first dispersion surface for the rectangular lattice with , in figure 7 possesses such a parabolic profile near . It is the horizontal spacing of that governs the special behaviour for normal incidence at , illustrated in figure 11 (see equation (26)). For the vertical period , we observe neutral propagation of the incident plane wave, similar to that observed for the single grating in figure 3(b), and we also see this for for the square lattice. Since the total displacement and incident fields are visually indistinguishable, we provide the real parts of the coefficients in figure 11(a) and the real part of the scattered field in figure 11(b). It is clear that for this choice of the wave does not “see” this specific lattice with , .

Neutrality for oblique incidence

For illustrative purposes we consider a square lattice of pins with aspect ratio before moving back to the rectangular lattice with . In figure 12, a plane wave is incident at with the Bloch parameter in the -direction set to .

Figure 12: A plane wave is incident at on a lattice of 2000 gratings with , for , . (a) Modulus of the coefficients for both Wiener-Hopf (dots) with , and Foldy (solid blue). The first 50 are shown. (b) Real part of total displacement field.

The real part of the total displacement field is shown in figure 12(b), and it appears that the direction of the incident field is virtually unchanged.

Figure 13: (a) Isofrequency contours for the first surface for square lattice with , . The ambient medium’s contour for is the dashed circle . (b) Isofrequency curves for the 2nd surface. The dashed line has been drawn to correspond with Bloch waves in the vertical direction for , . The Brillouin zone is also shown. (c) Close-up of isofrequency contours around . (d) Star-like features of the dispersion curves for at and collection of light circles for the square lattice (note the multiple intersections of the light lines at Dirac-like points).

Another notable feature is that the stripes of the incident plane wave are replaced by spots of positive and negative intensity; this pattern is due to the individual scatterers and is linked to the neutrality arising from the Dirac-like point. Figure 12(a) illustrates the moduli of the scattering coefficients obtained using both the Wiener-Hopf (dots) and Foldy (solid blue curve) methods. The agreement is very good, virtually indistinguishable for the first fifty gratings, and for both approaches the coefficients decay to a nonzero constant supporting the neutral propagation that we observe in part (b). One other interesting feature is that this propagation is not the only wave action, the appearance of “spots” to the left of the crystal indicates that some reflection action is also present.

In figure 13 we use isofrequency diagrams to aid our interpretation of this neutrality effect. In part (a) we show the isofrequency contours for the first surface for the square lattice, and in part (b) the second surface since is in the vicinity of Dirac-like cones for both. The contours surrounding the point in the Brillouin zone in both parts (a) and (b) indicate that is a Dirac-like point. The lowest value contour near in figure 13(a) is , and in part (b) for the second surface is . Part (c) shows a magnified picture of the neighbourhood of , for which the contour is more visible. We also illustrate an alternative derivation of the contours in figure 13(d), where they take the form of star-like features at the Dirac points. Figure 13(d) also displays the high degree of degeneracy of the light lines at Dirac-like points such as , where we observe multiple intersections. Note that figure 13(b) features the projections of light cone intersections onto the -plane in the form of lines for the square lattice defined by similar to equation (28) in section 3.2.1.

The dashed line intersects both the ambient medium’s contour and these platonic crystal isofrequency contours precisely at the point M. Recalling that the group velocity should be perpendicular to these slowness contours, and in the direction of increasing , we would expect to see refracted and reflected waves at an angle of as illustrated in figure 12(b).

Interfacial waves for the second dispersion surface

Now we change the aspect ratio to to obtain a rectangular lattice as in section 3.2.3, but keeping all other parameter values the same. In figure 14(a) we plot the isofrequency contours for the second dispersion surface, ranging from around 3.9 to 4.65.

Figure 14: (a) Isofrequency curves for the second surface for a rectangular lattice with , . The angle of incidence , the dashed line and the ambient medium’s contour are also shown. (b) Real part of the displacement field for , , for the rectangular lattice.

We also add the incident wave arrow for and the dashed line = 3.1113. The intersecting line at cuts the contour in multiple places meaning that it is more difficult to predict the behaviour of the system, although intersections corresponding to the group velocity being directed towards the interface from the crystal can be ruled out since the only incoming energy is from the incident medium [13].


Figure 15: (a) Solutions of the eigenvalue problem for a system of three gratings with and separation for . (b) Normalised transmitted energy versus for the triplet’s -1 diffraction order (solid black curve) and triplet’s 0 order (dashed red curve) for . (c) Normalised reflected energy versus for the triplet’s -1 diffraction order (solid black curve) and first pair’s 0 (dashed red curve) and -1 (thick blue) orders for . (d) Real part of displacement for , , .

The intersection at A indicates that the system supports transmission action in the form of refracted waves at an oblique angle less than . This transmission action is visible in figure 14(b), although at an apparently lower intensity than the reflection observed. Another interesting feature of the real part of the displacement field plotted in figure 14(b) is a trapped wave close to the edge of the half-plane system. This skew-symmetric mode is localised within the two channels formed by the first three gratings, which is an example of wave-guiding, with the reflection pattern incorporating the third grating as well as the first.

It appears that the triplet’s trapped mode is odd, so that in effect the central grating is not seen, and this demonstrates again the connection with finite-grating stacks discussed in section 3.2, although here we have coupling between the diffraction orders 0 and -1. We explain this coupling using figure 15, where we show the solutions to the eigenvalue problem for a set of three gratings with period , separation and .

Referring to figure 15(a), the first solution at is due to the zeroth order only, whereas the solution at arises for a mixture of 0 and -1 orders. This value of is very close to explaining the trapped mode we observe in figure 14(b), an effect we optimise by setting the frequency of the system with .

We show transmitted and reflected energy profiles for the corresponding transmission problem for the first two and three gratings in figures 15(b, c). In part (b), the solid black curve denotes the transmitted energy due to the -1 order for the triplet, explaining why the localised mode helps to support transmission through the rest of the half-plane system. The dashed red curve in part (b) represents the transmitted energy for the 0 order, emphasising that the -1 order dominates. This means that the first pair has to be able to facilitate transmission of the -1 order to reach the third grating. This is illustrated in figure 15(c) where reflected energy is plotted versus . The thick blue curve shows the reflection due to the pair’s -1 order is around meaning that there is at least transmission to feed the third grating, which in turn helps to feed the rest of the system.

Figures 15 (a-c) show that the maxima for the triplets arise for rather than . Therefore in figure 15(d) we plot the real part of the displacement field for the same rectangular lattice for , but with and the corresponding value of to support Bloch waves along the vertical gratings. The results are similar to those observed in figure 14(b) except that the transmission pattern within the pinned system is slightly more intense than that observed for .


Recalling the second surface’s isofrequency diagram 14(a) there are several notable features similar to those observed for the first surface. There is a parabolic profile along for and an inflexion at for where the contours change direction. Using the parabolic profile we are able to demonstrate an example of channelling in figure 16, where a plane wave with oblique angle of incidence is bent to travel parallel to the horizontal axis.

Figure 16: (a) Moduli (solid black) and real parts (dashed blue) of the coefficients for 4000 gratings with , , . The first 100 of 4000 gratings with , are shown. (b) Real part of the total displacement field.

The wave-vector diagram is used to identify isofrequency contours parallel to for the operating frequency for the rectangular lattice with , ; this ensures that the refracted waves will travel in a direction normal to the contours as illustrated by the point B in figure 14(a). The corresponding values for and are 2.2 and respectively. The resultant coefficient and displacement field plots are shown in figure 16.

Figure 16(a) shows very high values for both the real parts and the moduli of the coefficients; evidence that the system supports propagation. This is illustrated in part (b) where the reflection action is minimal, but the propagation through the pinned system is clearly visible, along with the platonic crystal’s wave-guiding effect.

Figure 17: (a) Moduli (solid black) and real parts (dashed blue) of the coefficients for 4000 gratings with , , . The first 30 of 4000 gratings with , are shown. (b) Real part of the total displacement field

A similar but more intense effect is observed for at with and . The resultant propagation parallel to is illustrative of the uni-directional localised modes associated with such parabolic profiles, and can be rotated through by swapping the horizontal and vertical periods of the lattice. Thus for a lattice of 4000 gratings with period and spacing , we show an interfacial wave in figure 17(b) propagating in the vertical direction parallel to , as predicted by the corresponding isofrequency contour, and very little wave action inside the pins as predicted by the decay of the coefficients in figure 17(a).

Higher frequencies - interfacial waves

It is clear that as the frequency increases, a more complicated picture emerges. The third surface illustrated in figure 6 and magnified here in figure 18(c) possesses many important features.

Figure 18: A plane wave is incident at on a lattice of 2000 gratings with , for , . (a) Moduli of the coefficients with the first 50 shown for both Foldy (solid blue curve) and Wiener-Hopf (dots). (b) Real part of the total displacement field. (c) Magnification of the third dispersion surface. (d) Real part of the total displacement field for , , .

At the point the contours change direction by for and there is a Dirac-like point at for . This narrow frequency window supports an array of interesting wave phenomena, similar to those discussed for the first surface in section 3.2.3.

We begin by considering a neighbouring frequency illustrated in figure 18(c) where we show a collection of isofrequency contours for the platonic crystal’s third dispersion surface, as well as the contour for the homogeneous biharmonic plate medium - an arc of the dashed circle defined by for the case . We also refer again to equation (28) to explain the appearance of the light cone projections on the isofrequency diagram.

We seek an interfacial wave so we choose such that is tangent to the contour for the third surface (point A in figure 18(c)). The vertical arrow predicts the propagation direction of the refracted waves for this and corresponding . In this way we predict that some refracted waves will propagate perpendicularly to this contour, i.e. in a vertical direction parallel to the interface. We plot the real part of the displacement field for these parameter settings in figure 18(b), along with the moduli for the coefficients in figure 18(a) generated by both the Foldy and Wiener-Hopf methods.

Figure 18(b) shows low propagation within the pinned structure, and relatively low reflection, compared with figure 15(b) for example. However there is clearly a travelling wave localised within the first few gratings, albeit with a relatively long wavelength. The exponential decay of the moduli of the scattering coefficients is illustrated in figure 18(a). In figure 18(d) we consider an example for an arbitrary value of to highlight the contrast in behaviour for various and . For and , we show the real part of the total displacement field in figure 18 (d). Although there is increased intensity at the interface, as one would expect, the dominant behaviour is propagation and in almost the same direction as the incident wave. This is predicted by the direction of the blue arrows at B and on the contour close to the incident wave’s arrow in figure 18 (c).

Figure 19: A plane wave is incident at on a lattice of 4000 gratings with , . (a, b) Real part of total displacement field for , . (c) Comparison of the moduli of coefficients .

Finally we consider the case , . At the point in figure 18(c), the contours for support propagation in a direction normal to those of , which support propagation parallel to the interface similar to A for for in figures 18(b,c). Thus there is a point of inflexion for and it occurs for . For , the refracted waves propagate through the pinned system rather than along its edge, which is what happens for . However at the Dirac-like point ; , the waves exhibit a mixture of edge localisation and propagation through the system. This is because we have the third, fourth and fifth dispersion surfaces meeting at this frequency, leading to a combination of behaviours. The third surface 18(c) predicts propagation along the edge for contours while the fourth and fifth surfaces predict propagation parallel to , into the pinned lattice.

We observe these properties in figure 19 where parts (a,b) show the real part of the total displacement field for and . The amplitudes are extraordinarily large for where the isofrequency contours are about to change direction to support interfacial waves rather than propagation through the pinned system (see in figure 18(c)). This transition has occurred in figure 19(a) for where the localised interface mode clearly dominates any action inside the pinned region. The comparison of the coefficients in figure 19(c) also demonstrates the decay of the coefficients to zero inside the crystal for . However close to a Dirac-like point supports both edge localisation and wave propagation as illustrated in figure 19(b, c), where a wave is seen to propagate through the pinned lattice with the wavelength of the envelope function matching that of the coefficients.

4 Concluding remarks - interfacial waves and dynamic localisation

The theory and examples presented in this paper have identified novel regimes, typical of flexural waves in structured Kirchhoff-Love plates, for the case of a semi-infinite structured array of rigid pins in an otherwise homogeneous plate. The interface waveguide modes have been identified and studied here.

For a half-plane occupied by periodically distributed rigid pins, we have demonstrated an interplay between the transmission/reflection properties at the interface and dispersion properties of Floquet-Bloch waves in an infinite doubly periodic constrained plate. Specifically, we have analysed regimes corresponding to frequencies and wave vectors that determine stationary points on the dispersion surfaces, as well as Dirac-like points. For such regimes we have demonstrated that the structure supports localised interfacial waves, amongst other dynamic effects.

The localisation was predicted by an analytical solution, and a formal connection has been shown here between the doubly quasi-periodic Green’s function for an infinite plane and the system of equations required to obtain the intensities of sources at the rigid pins, which occupy the half-plane.


The authors thank the EPSRC (UK) for their support through the Programme Grant EP/L024926/1.


Appendix A Factorizing

A critical technical detail is the factorization = with respectively analytic inside and outside the unit circle. Laurent’s theorem gives


where is a circle of radius slightly larger than the unit circle, and is a circle of radius slightly smaller than the unit circle (see figure 2); the desired factorization is


We note here that is of the form for and , with . It is also assumed that are chosen so that is well defined.

Referring to the factorization (A2),

is evaluated for for some argument on the unit circle, with for . Here we outline two possible choices of regularisation, both of which give good results. Defining with , we have


Then we may write


The alternative way is to define and then determine