## 4 Spatial Acoustic Scattering  
Figure 4 compares the 3D hemispherical scattered sound pressures of the aforementioned QRM with \(N=5\) slits and the equivalent RLC surface impedance, ZEQ, with \(N=5\)  
Figure 3: Magnitude and phase of the reflection coefficient using the Transfer Matrix Method (blue, solid) and using the RLC approximation (orange, dashed) relating two different slit admittances within the QRM along with fitting errors (green, dotted), as a function of frequency in Hz. slits as well. The scattering distributions have been obtained through different methods, i.e., (i) the Rayleigh-Sommerfeld (R-S) integral for near-field scattering information, (ii) Finite Element Method (FEM) in COMSOL Multiphysics 5.3 TM, and (iii) Finite-Difference Time-Domain (FDTD) simulations. Numerical simulations in COMSOL were computed by installing the surface at the centre of a spherical domain filled with air, surrounded by a concentric perfectly matched layer (PML) with a far-field boundary condition at the boundary of the air domain to satisfy [PERSON]'s radiation condition. A similar setup was used in the FDTD domain, but with first-order Engquist-Majda [40] absorbing boundary conditions at the boundary of a cubic air domain large enough to ignore any erroneous reflections in the obtained responses. The FDTD simulation was calibrated to have less than 1% numerical dispersion error up to 8 kHz with 10.5 points per wavelength (PPW) using the usual Courant-Friedrichs-Lewy (CFL) condition for the simplest 3D Cartesian scheme (\(CFL=1/\sqrt{3}\)), for which stability is ensured [21; 41]. The time step, \(T_{s}\), considered in the simulations can be deduced from the CFL condition as \(T_{s}=X_{s}\times CFL/c\), where \(X_{s}\) is the spatial sampling and \(c=343\) m/s is the speed of sound in air. Figure 4 shows the theoretical and numerical solutions of the surface scattered sound energy integrated over a radius distance of 1 m, so that all datasets represent a finite framework. In the FDTD model, virtual microphones were positioned at 1 m around the surface, which considering its edge dimensions (\(x,y\)) [35 cm, 35 cm], should be sufficient to correctly depict the scattered field at a frequency \(f=2\) kHz. FEM results were obtained following a similar approach where the integration was performed over a spherical near-to-far-field boundary condition of 1 m radius. Overall, the polar plots displayed in Figure 4a,b show some variations between the theory and FEM simulations, with normalised diffusion coefficient values varying from \(\delta_{n,0}=0.61\) to \(\delta_{n,0}=0.55\). These can be explained through the divergence of theoretical assumptions with respect to a numerical solving of the wave equation. More specifically, in a theoretical framework, the surface impedance is considered locally homogeneous for each slit, a fact that may not be entirely true in numerical terms due to the potential  
Figure 4: Distribution of sound scattered pressure levels at 2 kHz for a Quadratic Residue Metadditive (QRM) and the equivalent surface impedance RLC circuit (ZEQ) according to **(a)** Rayleigh-Sommerfeld (R-S) integral, **(b,d)** Finite Element Method (FEM), and **(c,e)** Finite-Difference Time-Domain (FDTD). \(\delta\) insets represent the diffusion coefficient of the scattered distributions. evanescent coupling between slits, a phenomenon not taken into account in the theory. As the scattered sound field is highly dependent on the distribution of the surface's reflection coefficient, slight variations in the polar distributions can thus be expected. Nonetheless, it can be seen that the global shapes of the QRM sound scattering distributions are sensibly similar one to another, which can be confirmed by their close autocorrelation values. In addition, Figure 4b represents quite well the expected main dip that defines the Quadratic Residue sequence at \(\theta\approx-20^{\circ}\) in the \(\phi=[-90^{\circ}:90^{\circ}]\) elevation axis. The main axial and lateral energy lobes are also quite well represented.