## 3 RLC Circuit Impedance Boundary Formulation  
The magnitude and phase of the surface impedance of a metadifuser are marked by the inherent resonances of the structure, i.e., if it has one or more resonators it should show one or more fundamental resonant peaks or phase shifts in its surface impedance. In the case of the QRM, a highly resonant structure with two HRs per slit is displayed, which lends itself to a boundary formulation based on a combination of second-order resonators. Example admittances for two slits of the QRM are shown in Figure 2. The two slits displayed in Figure 2 are both loaded with two HRs but with different slit and HR dimensions. This results in different surface impedances for each slit, thus presenting different resonance peak amplitudes and phase shifts. Figure 2: Normalised magnitude and phase of two example surface admittances (**left** and **right**) of two different slits within a QRM as a function of frequency in Hz, with identified resonances marked in grey dots. \(Y_{0}=1/\rho_{0}c_{0}\) represents the admittance of air. One useful passive formulation for a given slit impedance consists of a parallel set of series RLC circuits (after [28]), each consisting of one resistor (R), one inductor (L) and one capacitor (C), where all RLC parameters are non-negative and real-valued. The admittance of this structure in the Laplace domain is given by  
\[Y(s)=\sum_{b=1}^{B}\frac{1}{Z^{(b)}}\,,\quad Z^{(b)}=L^{(b)}s+R^{(b)}+1/(sC^{(b )}), \tag{4}\]  
where \(L^{(b)}\), \(R^{(b)}\), and \(C^{(b)}\) are non-negative constants for each branch. For our analyses and optimisations, we can limit \(s\) to \(s=j\omega\), and where \(B\geq 1\) is the number of RLC branches. The associated impedance of the RLC circuit is then simply \(Z(j\omega)=1/Y(j\omega)\). With this impedance boundary formulation we aim to fit the surface impedance of each slit of the QRM with an equivalent circuit made up of a set of resonances with non-negative RLC parameters. The passivity of this structure is then preserved in the discrete FDTD setting through the choice of the bilinear transform as a discretisation method [34]. This impedance model can be seen as an extension of various simpler frequency-dependent boundary models presented in the context of FDTD methods for room acoustics [35; 36; 37]. A methodology for fitting RLC parameters is described in [28], which consists of: (i) identifying resonances in an admittance by their peaks in the admittance magnitude, and (ii) estimating half-power bandwidths for each resonance from admittance magnitude. From those estimates RLC parameters may be identified for each resonance. This is followed by a global optimisation over the RLC triplet parameters using a Nelder-Mead optimisation. That approach works well when admittance peaks are well-separated, but in general peaks in admittance magnitude data can be difficult to identify, and furthermore half-power bandwidths can be hard to estimate from admittance magnitude data alone. This is especially true in the example admittances shown in Figure 2. Other ways for approximating the weight and resonant frequencies in viscoacoustic problems were also previously reported [38]. In this study, we use an approach which is based on making use of admittance phase information, and derivatives thereof, to identify resonance parameters. It can be observed from Figure 2 that peaks in admittance magnitude are linked to inversions in the phase response. More specifically, we know that the phase response of an individual-series RLC circuit admittance goes to zero at its resonant frequency and displays a negative slope at that frequency. Additionally, regarding the slope of the phase at the resonant frequency, one can derive from Equation (4)  
\[-\left.\frac{\partial\angle Y(i\omega)}{\partial\omega}\right|_{\omega=\omega _{0}}=\frac{2}{\Delta\omega}=\frac{2L}{R}, \tag{5}\]  
where \(\omega_{0}=1/\sqrt{LC}\).