## 3 Results  
As derived in the above section, relations (27) and (28) represent two phase unwrapping approaches based on the K-K relation (4) and the Hilbert transform (6), respectively. We named PUNW\({}_{KK}\) and PUNW\({}_{HT}\) the related numerical procedures. Regarding the PUNW\({}_{KK}\) method, the numerical evaluation of the involved K-K integral has been performed using the Simpson quadrature method described in [38], while for the PUNW\({}_{HT}\) procedure, the Hilbert Transform has been computed by using the FFT technique described in [34]. For both methods, the parameter \(p\) has been determined, forcing it to assume an integer value selected to ensure the continuity of \(arg_{\pi}(\cdot)\)[36]. To assess their numerical performances, we compared these with those provided by the numerical evaluation of the K-K relation (4) using the same numerical integration algorithm adopted for PUNW\({}_{KK}\)[38] and the phase unwrapping PUNWOS method [36] using the implementation described in [30]. To this end, we have considered recovering the effective refractive index \(n_{eff}(\omega)\) of two theoretical double negative MM slabs, in which relative permeability \(\mu_{r}(\omega)\) and permittivity \(\epsilon_{r}(\omega)\) are described by the following Lorentzian model [3]. \[\mu_{r}(\omega) =\mu_{\infty}+\sum_{n=1}^{np}\frac{(\mu_{s,n}-\mu_{\infty}) \omega_{0m,n}^{2}}{\omega_{0m,n}^{2}-\omega^{2}+i\gamma_{m,n}\omega} \tag{29}\] \[\epsilon_{r}(\omega) =\epsilon_{\infty}+\sum_{n=1}^{np}\frac{(\epsilon_{s,n}-\epsilon _{\infty})\omega_{0e,n}^{2}}{\omega_{0e,n}^{2}-\omega^{2}+i\gamma_{e,n}\omega} \tag{30}\]  
The parameters of the two models are reported in Table 1. This modeling choice is justified because: (_i_) Lorentzian models describe the dispersive nature of split-ring resonators and strip wires [3], and (_ii_) the superposition of Lorentzian functions describes every causal function [39]. The MATLAB environment has been employed to code all methods. All numerical experiments were performed on a 2.0 GHz Intel Xeon workstation. The percentage error (PE) is defined as follows:  
\[\frac{\|n_{ex}(\omega)-n_{eff}(\omega)\|_{2}}{\|n_{ex}(\omega)\|_{2}}\% \tag{31}\]  
where \(\|\cdot\|_{2}\) is the 2-norm [40], \(n_{ex}(\omega)=\text{Re}[\sqrt{\epsilon_{r}(\omega)n_{r}(\omega)}]\) and \(n_{eff}(\omega)\) are the exact and recovered refractive indexes, respectively, which have been used as criteria to evaluate the efficiency of the considered recovering procedures. Moreover, computational time t has been evaluated. However, this last datum has to be considered as only indicative of the order of magnitude rather than the exact time involved in each considered retrieval procedure since such knowledge would require the development of numerical codes optimized for a given CPU in order to achieve the best possible computational performances [41]. Figure 1 shows path \(\gamma(\omega)\) traced by Equation (8) for the first slab, with \(d_{eff}=180\) nm, which is characterized by a single pole Lorentzian model for pair \((\epsilon_{r}(\omega),\mu_{r}(\omega))\). As pointed out in [30], this trend for \(\gamma(\omega)\) confirms that this is an electrically thick slab because \(\gamma(\omega)\) repeatedly intersects ray \(\mathbb{R}_{\pi}\), i.e., the branch cut for the principal complex logarithm \(log(\cdot)\)[30; 32]. Accordingly, the complex exponential function in equation \(e^{(\cdot)}\) in (1) lacks univalence for the considered case, thus indicating that \(n_{eff}(\omega)\) must be evaluated via the right-inverse operation [30]. Table 2 shows the values of PE and t as a function of the number of sample points \(n_{s}\) in frequency band \([0,1]\) PHz for all the considered recovering methods.