## 2 Governing Equations and Methods
### Transient GPU-Accelerated FE Models  
The underlying finite element model's spatial domain is illustrated in Figure 3a for the case of an incident Gaussian beam. Various dimensions of the metamaterial slab were considered but the case of ten rows of sixty identical unit cells is shown in Figure 3a, covering horizontal and vertical distances of 60 mm and 10 mm, respectively. Absorbing layers, to reduce the effect of unwanted reflections from the sides of the truncated two-dimensional homogeneous part of the domain, are applied to the horizontal and vertical extremities of the domain. The thickness of the absorbing layers is chosen to be \(>3\lambda_{\mathrm{p}}\), following the recommendations of previous literature [33] to ensure that unwanted reflections are minimised; the notation \(\lambda_{\mathrm{p}}\) denotes the longitudinal, or compression, wavelength. The dependence of the absorbing layer parameter value, on the longer wavelength scale, accounts for the mode conversions that occur as the incident shear wave interacts with the metamaterial and its constituent subunits. The wavelength is determined by the incident centre frequency, \(f=\omega/2\pi\), and the material parameters of the media. Aluminium material parameters (Young's modulus, Poisson's ratio, and density) are assumed here: \(E=70\) GPa, \(\
u=0.33,\rho=2700\) kg/m\({}^{3}\) to define compression and shear wave speeds of \(c_{\mathrm{p}}=6198\) and \(c_{\mathrm{s}}=3122\) m/s. The two homogeneous parts of the domain in Figure 3a, above and below the metamaterial slab, are defined by these aluminium material parameters. The central portion was also initially defined in the same way but the cross-shaped holes are excised to create the periodic array of unit cells. Zero traction boundary conditions (see Equation (3)) are applied to the edges of the crosses, whose dimensions are illustrated in Figure 1a and Table 1.  
Each unit cell is a 1 mm \(\times\) 1 mm square that contains a symmetrical cross (resembling the Swiss cross) with length and width 0.9 mm. Therefore, the spacing between the left and right ends of neighbouring crosses is 0.1 mm. As shown in Figure 1a, all sublengths of this study's fundamental cross are equal, i.e., 0.3 mm, whereas the Swiss cross (also known as the emblem of the Red Cross) possesses a sublength ratio of 7:6, for which the central square has relative length 6. The initial studies investigated for this article were periodic arrays of these regular crosses, with slab sizes varying from \(10\times 60\) slabs to \(24\times 60\), for which the first number represents the number of repeating rows. The domain was meshed with element length \(\Delta x=0.02\) mm using the explicit time-domain finite element software package Pogo [24], and its in-built meshing tool, pogoMesh. The element size is minimised to improve the  
Figure 3: (a) Finite element model domain for shear wave Gaussian beam incidence. The angle \(\alpha\) is \(35^{\circ}\) and the slab of metamaterial, sandwiched by homogeneous aluminium media, consists of 600 cross-shaped holes. Absorbing layers of length 8 mm are applied to all sides of the domain. (b) Alternative point-source excitation for shear waves. The weighting applied to each of the six excitation nodes is highlighted by the directions of the arrows. accuracy of the results observed. It is well known that for two-dimensional explicit time domain FE studies [24; 34], thirty elements per wavelength is more than sufficient for convergence of the results. Since a range of frequencies were used here, the mesh size was chosen to ensure convergence for all cases considered. The highest centre frequency in the examples that follow is \(3.3\,\mathrm{MHz}\), which corresponds to a shear wavelength of \(\lambda_{\mathrm{s}}=0.95\,\mathrm{mm}\), equating to close to fifty elements per shear wavelength, ensuring convergence for all examples presented.