## 2 Governing Equations and Methods
### Transient GPU-Accelerated FE Models
Figure 4a illustrates the mesh of plain-strain, linear, 3-noded triangular elements around a cross-shaped hole; the pogoMesh function generates an irregular mesh in the immediate vicinity of the edges and vertices, but reverts to a regular form as the distance from the interface increases. The size of the FE model varies with the slab of metamaterial investigated. For example, a \(10\times 60\) example, with beam angle \(\alpha=35^{\circ}\) varies from \(-26\) mm to \(36\) mm in the vertical direction, and \(-40\) mm to \(40\) mm in the horizontal dimension, producing a domain consisting of around \(25\) million triangular elements and degrees of freedom (DOFs) for the two-dimensional model. A larger domain is required to accommodate a smaller beam angle, or a point source at distance. For example, the case of a point source located \(40\) mm from the top of the slab requires a domain with more than \(50\) million DOFs, ranging from \(-40\) mm to \(40\) mm in the horizontal direction, and from \(-56\) mm to \(56\) mm in the vertical direction. The use of Pogo [24], a GPU-driven explicit transient FE software package, has enabled a step change in modelling capability for the elastic metamaterial considered here. Standard CPU-driven software packages for time domain computations begin to show instability for model-sizes \(>10\) million DOFs, and simulations become increasingly inefficient with respect to run-times. In contrast, several examples considered here exceeded \(50\) million DOFs with no stability issues and run-times of less than ten minutes on a standard PC with a single Nvidia GTX 1080 Ti GPU card possessing \(11\) GB of memory. The results are also highly accurate for very large and complex models; the three-dimensional studies of elastic wave attenuation and scattering in polycrystalline solids [26; 35] investigate models with up to \(1\) billion DOFs. Models that would take several days for standard CPU-driven transient schemes can be run in a matter of hours using parallel GPU-implementation. The underlying excitation is a \(200\)-cycle Hann-windowed tone-burst (numerical investigations were initially conducted to optimise the number of cycles) with centre frequencies selected from the higher frequency dispersion surfaces in Figure 2. Two types of shear wave excitation were considered, a Gaussian beam or point source. The Gaussian beam used a \(5\) mm source line, highlighted with red markers in Figure 3, wherein the underlying signal was applied to each source node with appropriate weighting to steer the beam at the chosen incident angle \(\alpha\) (see Figure 3a). The point-source excitation method is illustrated in Figure 3b. Similarly, appropriate weightings were applied to the relevant nodes to produce an incident shear wave, whose homogeneous medium wavelength is determined by \(\lambda=c_{\mathrm{s}}/f\), where \(f\) denotes the centre frequency of the excitation. Note that a phase difference of \(\pi\) is applied to nodes diametrically opposite to one another (indicated by arrows of the same colour) in Figure 3b. In order to generate a compression wave point source, the arrows in Figure 3b would be rotated by \(\pi/2\) in the anticlockwise direction, to obtain the required weightings. Each simulation was run for \(T=120\)\(\mu\)s, with the excitation signal length defined by the ratio of number of cycles to the frequency. For explicit time-domain numerical methods, it is essential that the Courant-Friedrichs-Lewy convergence condition [36] is satisfied. The finite element simulation time step \(\Delta t\) was determined by the choice of Courant number \(C\) and the equation [36]  
\[\Delta t=C\Delta x/c_{\mathrm{p}}. \tag{7}\]  
The Courant number \(C\) must be less than one, and it is recommended to be lower than \(0.7\), so our choice of \(C=0.3\) ensures convergence for accurate results. The mode-converted compression wave speed, rather than shear wave speed, is used in Equation (7) to ensure that the time steps are sufficiently small to capture the wave propagation with optimal accuracy.