Analytic Continuation, Phase Unwrapping, and Retrieval of the Refractive Index of Metamaterials from S-Parameters

The heuristic homogenization approach is intensively employed to characterize electromagnetic metamaterials (MMs). The effective parameters are extracted within this framework using the Nicolson–Ross–Weir (NRW) method. Special attention must be devoted to handling this procedure because of the branch ambiguity issue affecting it, i.e., the lack of uniqueness in the evaluation of the effective refractive index neff rooted in the use of the multivalued complex logarithm to invert the Airy–Fresnel relation. Over the years, several techniques based on the phase-unwrapping approach have been introduced, but without any theoretical justification. In this paper, we aim to clarify the theoretical connection between the phase unwrapping method and the analytic continuation theory framework. Furthermore, three-phase-unwrapping approaches, which descend directly from the theory we discussed, are compared to identify which approach is best suited to reconstruct the complex refractive index of metamaterials when the NRW method is applicable.


Introduction
Retrieving the electromagnetic parameters of materials is a task of primary importance in many research fields, such as microwave engineering, electromagnetic compatibility, and bioelectromagnetics, to name a few [1][2][3].In the last decade, researchers have shown great interest in developing particular three-dimensional artificial materials, usually made up of a lattice of metallic resonant inclusions arranged in a dielectric host medium called metamaterials (MMs), given obtaining devices endowed with exceptional operating performances [4].To characterize MMs, the so-called heuristic homogenization approach, developed for the first time in [5,6], is the procedure most commonly employed by practitioners and researchers in the field [7][8][9][10][11][12][13][14].In this peculiar framework, it is assumed that an MM, within its operating frequency range, can be considered analogous to a homogeneous medium, called the effective medium, for which the electromagnetic properties are described by an effective electric permittivity ϵ e f f and magnetic permeability µ e f f [5,6].Such effective permittivities are extracted by processing the numerical or measured scattering parameters of the MM at hand by using the Nicolson-Ross-Weir (NRW) method [5].Despite the intrinsic plainness of the NRW method, this recovering procedure often results in nonlocal ϵ e f f and µ e f f [15,16], and this poses severe problems regarding its scope of applicability [15,16].Nonetheless, the complex refractive index N e f f provided by the NRW approach can still be used for recovering the effective parameters for one particular class of MMs, called Bloch lattices [15][16][17].For these MMs, the effective parameters recovered by using N e f f in conjunction with the concept of the Bloch impedance Z B are consistent, obeying the locality constraints [17].More precisely, for this class of structures, the description in terms of ϵ e f f and µ e f f can be extended far beyond the quasi-static limit, reaching the frequency range in which the MM operates [17].To evaluate N e f f correctly, special attention has to be devoted to its recovering to provide correct results because of the branch ambiguity issue affecting the NRW method, i.e., the lack of uniqueness of Re[N e f f ], and the effective refractive index n e f f , which is rooted in the use of the multivalued complex logarithm LOG(•) to invert the Airy-Fresnel relation [18].In the literature, to face this ambiguity problem, several strategies grounded on the phase-unwrapping approach have been introduced [19][20][21], although without providing any theoretical justification for their use [22].Only recently, the link between phase unwrapping and the operation of analytic continuation of holomorphic functions has been demonstrated [22].Following a previous couple of works [22,23], in this study, we aim to better clarify the connection between the phase unwrapping method and the analytic continuation framework and identify which way is best suited to reconstruct the refractive index of metamaterial structures.The paper is organized as follows: Section 2 provides the elements of theory needed to relate analytic continuation and phase unwrapping.In particular, in Section 2.1, we summarize the NRW approach; in Section 2.2, we discuss the link between phase unwrapping and the analytic continuation theory and how the branch issue problem affecting the NRW method can be overcome using the phase-unwrapping approach.In Section 2.3, we connect the Kramers-Kronig integrals to phase unwrapping using the concept of the Riemann surface.In Section 3, the numerical performances of the phase-unwrapping procedures introduced in the above subsections are investigated.To this end, three different samples of a µ-negative MM composite realized via a three-dimensional array of LiTaO 3 spheres embedded in the free space have been considered.A numerical comparison of the phase unwrapping findings with those provided by the K-K relations and the Maxwell-Garnett theory is performed to rank the methods.Finally, in Section 4, conclusions are drawn.

The NRW Approach for Recovering the Complex Refractive Index
The Nicolson-Ross-Weir (NRW) method is a standard technique for recovering the permittivity and permeability of linear, isotropic, and homogeneous electromagnetic media [24] through closed-form relations derived considering the reflection-transmission phenomenon involving a plane wave normally incident on a slab of material with thickness d, placed in free space.Due to its inherent conceptual simplicity, the NRW method has also been applied as a homogenization method for metamaterials, as discussed extensively in [5,6].According to this approach, the complex effective refraction index N e f f of an MM can be recovered by solving the following relation (the time-dependence convention e −iωt is adopted) [22,24]: where k 0 , S 11 , and S 21 are the the free space wave number, and the S-parameters of the MM under study.The reflection coefficient between the effective medium slab and the free space, R, appearing in relation ( 1) is as follows: where z represents the MM effective intrinsic impedance.Using the complex logarithm LOG(•) = log(•) + 2pπi to solve (1), in which the term log(•) denotes the principal logarithm [25], and the integer p ∈ Z represents the branch index, we obtain the complex refractive index N e f f as follows: From (3), the expressions of the effective extinction factor κ e f f , the effective refractive index, and n e f f read as follows: )

Analyticity and Phase Unwrapping
To avoid the lack of unicity of the refractive index n e f f , when solving (1), it is necessary to take into account some preliminary considerations about the invertibility of the complex exponential function e (•) .More precisely, relation (1) will be invertible, and its inverse will be unique if and only if e (•) is univalent [22,25].Unfortunately, this property is related to the term iN e f f k 0 d, which is the unknown of the problem at hand.Accordingly, we must treat e (•) as not univalent, and since the inverse function does not exist in this case, we have to solve (1), computing its right inverse [22].The following theorem provides us with a method for its evaluation [22]: fulfilling the constraint condition Im(L(γ(a))) = θ 0 with θ 0 ∈ R and e iθ 0 = γ(a) |γ(a)| .
The detailed proof of the above theorem is given in [22].In the following, we will summarize the rationale and highlight the most essential points.First, the solution of the equation is calculated by superimposing the path γ(t) with a series of disks, thus realizing the so-called path covering [25].Each disk is the domain of an appropriate analytic logarithm.In this way, the analytic continuation operation is set up.Second, the analytic logarithm solution of ( 6) results as the appropriate superposition of the analytic logarithms defined within each disk, for which the phase term is made appropriately continuous as it passes from one disk to another.The operation essentially entails the restoration of the continuity of the imaginary part of L(•) that, as an important consequence of Theorem 1, can be obtained through a phase-unwrapping approach.In the particular case of Equation ( 1), when e (•) is right-inverted using the above result, it is fundamental to use a suitable constraint condition to obtain a unique significant solution from a physical point of view.Considering that γ(t) ≜ S 21 /(1 − S 11 R) with t = ω, the constraint Im(L(γ(a))) = θ 0 has to be specialized as follows: which ensures the causality of the refractive index [22].In fact, this choice ensures that Im[L(γ(t))] is an odd function.In this way, we have that L(•) results as the Fourier transform of a time-domain causal physical quantity [22], and that the principal logarithm log(•) can be used as a first element of the chain of analytic function elements realizing L(•).Having calculated the right inverse L(•), N e f f is given by the following: The computational results provided by Theorem 1, summarized in Algorithm 1, justify rigorously the use of this type of approach in the literature, usually implemented by using the Oppenheim and Schafer algorithm exploited in the signal processing field [26] (see Appendix A for more details on this point), for recovering the MM's electromagnetic parameters [22].An algorithm derived from the phase-unwrapping procedure of Oppenheim and Schafer that, in principle, can locate and identify the crossings between γ(t) and the principal logarithm branch cut in a more precise way has been developed in [27].Its pseudo-code is reported in Algorithm 2. A sampled version of γ(t) ≜ S 21 /(1 − S 11 R), [γ(t)] is given as its input.The algorithm checks if [γ(t)] crosses the branch cut of log(•), comparing the position of the points [γ(t)] j−1 and [γ(t)] j in the complex plane C − 0. The algorithm updates the branch index p accordingly, and computes the numerical value of [L(γ(t) j )] = log([γ(t)] j ) + i2π p at point [γ(t)] j .If an ambiguity arises, the algorithm stops, the path is re-sampled, using N = 2 np+1 points to this end, and the procedure is performed again.end if 23: end for

Riemann Surfaces and Phase Unwrapping
An alternate and more commonly employed method to overcome the ambiguity afflicting n e f f is to resort to the Kramers-Kronig relations [28][29][30].Because in any physical system the cause cannot precede the effect, the effective permittivities ϵ e f f and µ e f f must obey the causality principle [31,32].However, the complex refractive index is defined as follows: in which at a first glance it could seem that N e f f lacks analyticity, because of the presence of branch points in the upper half-plane (UHP) of C due to the zeros of the terms ϵ e f f and µ e f f in this region.Despite this, as demonstrated in [31,32], the term (9) does not have any branch point in the UHP in the case of passive media.This fact allows relating the real and the imaginary parts of N e f f (ω) via K-K relations and since the term κ e f f (ω) is uniquely determined from relation (4), it allows to determine n e f f (ω) without ambiguity as follows: To understand how the relation ( 9) can be exploited in a phase-unwrapping scheme, we have to reconsider the NWR Equation ( 1) from a more abstract point of view.As discussed in [23], the set of all the analytic logarithms L(•) of ( 6) makes up a special mathematical object called a global analytic logarithm L. A particular domain of definition is related to L: the Riemann surface S(L) [23].Its structure is composed of a numerable set of replies of Ċ, S p (L), overlayed on top of each other and sorted by the ascending index p, suspended above Ċ, and glued to each other along the slit extending from 0 to −∞, which each sheet owns [23], as shown in Figure 1.If z ′ is a point of Ċ and z is the point belonging to S q (L), for the q-th sheet of S(L), which lies above the point z ′ of Ċ, the value taken by L on z, L z′ , is given by [23]: Based on the above, if γ is the copy on S(L) of the path γ(t) lying in Ċ, we have from (11) that between the values taken by L(•) at z ′ ∈ γ(t) and the values assumed by L at z ∈ γ, the following relation must hold [23]: where, for each z ∈ γ, the value assigned to the integer parameter p is the index q of the Riemann sheet S q (L) on which z lies.Relation (12) can be used to compute L(•) on γ(t) ≜ S 21 /(1 − S 11 R).Substituting ( 8) into (12), we have the following: The unknown indexes q that localize the Riemann sheets S q (L) where the term S 21 /(1 − S 11 R) lies can be evaluated considering that from (8) we have the following: Inserting ( 10) into ( 14), we obtain for the q values as follows: However, from a computational point of view, considering the unavoidable numerical errors that are carried out in the evaluation of the K-K integral (10), the evaluation of the q values is better accomplished through the following minimization problem [26]: ) in which the q integers are those minimizing the error, ϵ err , in the absolute value of the difference between the argument of the principal logarithm and the term n e f f k 0 d.The two Equations ( 16) and ( 17) can be regarded as a numerical integration-like phase-unwrapping approach [26].As a matter of fact, the above q values can be recognized as the integers that allow the phase term arg π (S 21 /1 − S 11 R) to be continuous, i.e., to unwrap the argument of the principal logarithm, and thus provide the determination of the argument of L(S 21 /1 − S 11 R), the imaginary part of the right-inverse of which we are searching.Algorithm 3 reports its pseudo-code.

Numerical Results
To compare the performance of the phase-unwrapping approaches described in Sections 2.2 and 2.3 (denoted as Alg1 (phase unwrapping), Alg2 (plane phase unwrapping), and Alg3 (numerical integration-like phase unwrapping), respectively), we have considered the recovering of n e f f for a composite realized by a three-dimensional array of lithium tantalate (LiTaO 3 ) spheres, considered infinitely large along the x − y axis, and with finite thickness along z, embedded in free space [33].This is a µ-negative MM for which effective parameters can be analytically checked using the Maxwell-Garnett rules from mixtures [33].The geometry of the array is reported in Figure 2. Three structures, characterized by a different thickness d (which is an integer multiple of the cell size s; see Figure 2), were numerically simulated using the software Ansys HFSS, running on an Intel Xeon DP E5405 Quad Core-based workstation.The radius of the LiTaO 3 spheres was r 0 = 4 µm, and the unit cell size was s = 10 µm [33].The phase-unwrapping methods were programmed using MATLAB R2023b.The computation of the K-K integral (9) was implemented using the Euler-Maclaurin method described in [34] (the method is denoted as KK(B) in all figures and tables).Figure 3 shows the results for the case d = 30 µm.The rapid change of the plot of the phase of the S-parameters sampled with n s = 512 equispaced sampling points in the band (0, 5) THz, shown in the top left part of Figure 3, suggests that n e f f is affected by the branch ambiguity.This consideration is confirmed by the behavior of arg π (•), reported in the top left part of the same figure, which is clearly discontinuous, and needs to be unwrapped to compute the above parameter.Furthermore, the top right part of Figure 3 shows the comparison between the refraction index n e f f computed (i) using the Maxwell-Garnett theory [33], (ii) using the K-K integral, and (iii) using the phase-unwrapping procedures introduced above.We point out that all the considered phase-unwrapping procedures agree with each other and provide results that are close to the n e f f obtained by applying the Maxwell-Garnett theory.By contrast, the Kroenig-Kramers integral method reconstructed a refractive index profile that is not comparable with the result provided by the unwrapping methods (see the bottom left part of Figure 3 for details on this point).In the bottom-right part of Figure 3, the values assumed by the index q are reported as a function of the frequency.From these data, we can see that the Riemann sheets involved in the unwrapping process are S 0 (L) and S 1 (L).Table 1 reports the CPU time needed for each method for the considered case.The results show no particular superiority in terms of this parameter of one method over another.In Figure 4 are reported the results for the case d = 70 µm.As this structure is thicker than the previous one, we expect the term S 21 /1 − S 11 R to cross the branching line of the principal logarithm a more significant number of times than in the previous case.This fact is confirmed by the behavior of the phase of the S parameters, in this case sampled with n s = 1024 equispaced sampling points in the band (0, 5) THz.In fact, we can observe that it varies more rapidly than in the previous case.As a consequence, the term arg π (•) results to be more wrapped with respect to the case of d = 30 µm (see the top-left part of Figure 4.The unwrapping procedures provide results that agree with each other and with those provided for the previously considered less-thick structure (in fact, the refractive index of an assigned media is independent of its thickness).Also, in this case, the result provided by the K-K integral is less accurate (characterized by an oscillatory behavior) than that of the phase-unwrapping approaches provided, as clearly shown in Figure 4, in its bottom-left part, where the magnification of all curves in the (2, 3) THz frequency range is reported.S 0 (L), S 1 (L), and S 2 (L) are the Riemann sheets involved.The data reported in Table 2 show that the performances of the CPU time are comparable among methods.Finally, we considered an even thicker structure, with d = 130 µm.The number of sampling points used in this case was n s = 8046.The behavior of the argument of the S-parameters and of the argument of the principal logarithm, arg π (•), resulted in a more wrapped result than in the previous cases with d = 30 µm and d = 70 µm, as expected (see the top-left and the top-right parts of Figure 5, respectively).For this latter case, the results given by the phase-unwrapping procedures Alg1 and Alg2 are still in good agreement between them and close to the n e f f calculated using the Maxwell-Garnett mixing formulas.Regarding the Alg3 procedure, it can be noticed that the provided result is characterized by a discontinuity around f = 4.8 THz, a numerical error probably due to the fact that the KK integral is, for the considered case, affected by a large ripple, as shown in the bottom right of Figure 5.Nevertheless, the correction effect offered by the phase-unwrapping algorithm Alg3 on the poor result obtained by calculating the Kramers-Kronig integral is remarkable, as shown in the bottom-left part of Figure 5. S 0 (L), S 1 (L), S 2 (L), and S 3 (L) are the Riemann sheets involved in this last case.In terms of CPU time, the Alg3 results are more beneficial compared to the other methods.This result in terms of CPU time can certainly be attributed to the characteristics of this particular method, described in Section 2.2, and is now evident because of the more significant number of n s points used compared to previous cases (see Table 3).

Conclusions
Based on previous works [22,23], in this paper, after clarifying the connection between the phase unwrapping method and the analytic continuation theory, we have investigated the performances of different phase-unwrapping procedures.To this aim, we have conducted numerical experiments considering a µ-negative MM composite realized with a three-dimensional array of LiTaO 3 spheres, an artificial composite, intensively considered in the MM literature, for which effective parameters can be analytically checked using the Maxwell-Garnett theory.Numerical results show that the methods based on the simple unwrapping of the phase of the principal logarithm are adequate for all the considered cases, where the unwrapping method based on using the K-K integral can fail if the MM at hand is quite thick.However, this last method has the characteristic of correcting the numerical n e f f provided by the K-K integral relation (10), and this is a remarkable result.Regarding CPU time, all unwrapping methods are competitive, although the Alg3 method seems more advantageous for thicker structures.To conclude, we point out that the phase-unwrapping methods, already introduced in the MM literature [19][20][21] but without any theoretical justification, are well grounded in the analytic continuation theory and can be used without them being considered simply empirical techniques for the MM characterization.

Figure 1 .
Figure 1.A portion of the Riemann surface S(L) suspended on the unitary disk D ∈ Ċ.

Figure 2 .
Figure 2.An array of LiTaO 3 spheres embedded in free space: (a) sketch of the considered geometry; (b) unit cell containing a single LiTaO 3 sphere (s, cell size; r 0 , sphere radius).

Figure 3 .
Figure 3. Arrays of LiTaO 3 spheres, d = 30 µm; top-left: phase of the S-parameters; top-right: plot of n e f f ; bottom-left: magnification of n e f f over the (2, 3) THz band; bottom-right: Riemann sheets' q-values.

Figure 4 .
Figure 4. Arrays of LiTaO 3 spheres, d = 70 µm; top-left: phase of the S-parameters; top-right: plot of n e f f ; bottom-left: magnification of n e f f over the (2, 3) THz band; bottom-right: Riemann sheets' q-values.