Open AccessArticle

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

Department of Information Engineering, Infrastructures and Sustainable Energy, Mediterranea University, 89124 Reggio Calabria, Italy

Department of Civil, Energetic, Environmental and Material Engineering, Mediterranea University, 89124 Reggio Calabria, Italy

Dipartimento di Ingegneria Industriale e dell’Informazione, University of Pavia, 27100 Pavia, Italy

Author to whom correspondence should be addressed.

Sensors 2024, 24(3), 912; https://doi.org/10.3390/s24030912

Submission received: 6 January 2024 / Revised: 26 January 2024 / Accepted: 29 January 2024 / Published: 30 January 2024

(This article belongs to the Special Issue Innovative Devices and MEMS for Sensing Applications)

Download

Browse Figures

Versions Notes

Abstract

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

n_{e f f}

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.

Keywords:

phase unwrapping; scattering parameters; metamaterials

1. 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

L O G (\cdot)

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₃ 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.

2. Theory

2.1. 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]:

\begin{matrix} e^{i N_{e f f} k_{0} d} = \frac{S_{21}}{1 - S_{11} R}, \end{matrix}

(1)

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:

R = \frac{z - 1}{z + 1} w i t h z = \pm \sqrt{\frac{{(1 + S_{11})}^{2} - S_{21}^{2}}{{(1 - S_{11})}^{2} - S_{21}^{2}}},

(2)

where

z

represents the MM effective intrinsic impedance. Using the complex logarithm

L O G (\cdot) = l o g (\cdot) + 2 p π i

to solve (1), in which the term

l o g (\cdot)

denotes the principal logarithm [25], and the integer

p \in Z

represents the branch index, we obtain the complex refractive index

N_{e f f}

as follows:

\begin{matrix} N_{e f f} (ω) & = \frac{i}{k_{0} d} [l o g (\frac{S_{21}}{1 - S_{11} R}) + 2 p π i] . \end{matrix}

(3)

From (3), the expressions of the effective extinction factor

κ_{e f f}

, the effective refractive index, and

n_{e f f}

read as follows:

\begin{matrix} κ_{e f f} = - \frac{1}{k_{0} d} Re [l o g (\frac{S_{21}}{1 - S_{11} R})] \end{matrix}

(4)

\begin{matrix} n_{e f f} = \frac{1}{k_{0} d} Im [l o g (\frac{S_{21}}{1 - S_{11} R}) + 2 p π] . \end{matrix}

(5)

2.2. 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^{(\cdot)}

. More precisely, relation (1) will be invertible, and its inverse will be unique if and only if

e^{(\cdot)}

is univalent [22,25]. Unfortunately, this property is related to the term

i N_{e f f} k_{0} d

, which is the unknown of the problem at hand. Accordingly, we must treat

e^{(\cdot)}

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]:

Theorem 1.

Let

γ : [a, b] \to \dot{C}

be a continuous path. There exists a unique analytic logarithm

L (\cdot)

such that:

e^{L (z)} = z \forall z \in γ (t) t \in [a, b],

(6)

fulfilling the constraint condition

Im (L (γ (a))) = θ_{0}

with

θ_{0} \in R

and

e^{i θ_{0}} = \frac{γ (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 (\cdot)

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^{(\cdot)}

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:

\begin{matrix} Im {[L (\frac{S_{21}}{1 - S_{11} R})]|}_{a = ω_{0} = 0} = 0, \end{matrix}

(7)

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 (\cdot)

results as the Fourier transform of a time-domain causal physical quantity [22], and that the principal logarithm

l o g (\cdot)

can be used as a first element of the chain of analytic function elements realizing

L (\cdot)

. Having calculated the right inverse

L (\cdot)

N_{e f f}

is given by the following:

N_{e f f} = \frac{1}{i k_{0} d} L (\frac{S_{21}}{1 - S_{11} R}) .

(8)

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

l o g (\cdot)

, 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})] = l o g ({[γ (t)]}_{j}) + i 2 π p

at point

{[γ (t)]}_{j}

. If an ambiguity arises, the algorithm stops, the path is re-sampled, using

N = 2^{n p + 1}

points to this end, and the procedure is performed again.

Algorithm 1 Phase unwrapping

a : = 0

;

a r g_{α_{0}} (\cdot) : = a r g_{π} (\cdot)

;

γ_{t_{0}} : = \frac{S_{21} (a)}{1 - S_{11} (a) R (a)}

;

{[L (γ_{t})]}_{k = 0} = l n | γ (a) | + i a r g_{α_{0}} (γ (a))

;

5: for k:=1 to n do;

t_{k} : = ω_{k}

;

γ_{t_{k}} : = \frac{S_{21} (t_{k})}{1 - S_{11} (t_{k}) R (t_{k})}

;

p_{k} : = (a r g_{α_{k - 1}} (γ_{t_{k}})) - a r g_{α_{k}} (γ_{t_{k}})) / 2 π

;

a r g_{α_{k}} (\cdot) : = a r g_{α_{k}} (\cdot) + 2 π p_{k}

;

10:

[L (γ_{t}))]_{k} : = l n | γ_{t_{k}}) | + i a r g_{α_{k}} (γ_{t_{k}}))

;

11: end for

Algorithm 2 Plane phase unwrapping

a : = 0

;

{[γ_{t}]}_{0} : = \frac{S_{21} (a)}{1 - S_{11} (a) R (a)}

;

{[u_{t}]}_{0} : = Re ({[γ_{t}]}_{0})

;

{[v_{t}]}_{0} : = Im ({[γ_{t}]}_{0})

;

5: for k:=1 to n do;

t_{k} : = ω_{k}

;

{[γ_{t}]}_{k} : = \frac{S_{21} (t_{k})}{1 - S_{11} (t_{k}) R (t_{k})}

;

{[u_{t}]}_{k} : = Re ({[γ_{t}]}_{k})

;

{[v_{t}]}_{k} : = Im ({[γ_{t}]}_{k})

;

10: end for

11:

p : = 0

;

12: s_s

: =

¬ false;

13: for k:=1 to n do;

14:

{[L [γ_{t}]]}_{k} : = l o g ({[γ_{t}]}_{k - 1}) + i 2 π p

;

15: if

{[v_{t}]}_{k - 1} {[v_{t}]}_{k} < 0

then

16: if

{[u_{t}]}_{k - 1} \leq 0 \land {[u_{t}]}_{k} \leq 0

then

17:

p = p - sgn ({[v_{t}]}_{k - 1})

18: else if

({[u_{t}]}_{k - 1} \leq 0 \land {[u_{t}]}_{k} \geq 0) \lor ({[u_{t}]}_{k - 1} \geq 0 \land {[u_{t}]}_{k} \leq 0)

then

19: s_s

: =

¬ true;

20: Stop and Re-execute;

21: end if

22: end if

23: end for

2.3. 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:

N_{e f f} = \sqrt{ϵ_{e f f} μ_{e f f}},

(9)

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:

\begin{matrix} n_{e f f} (ω) - 1 = \frac{2}{π} P \int_{0}^{+ \infty} \frac{ω^{'} κ_{e f f} (ω^{'})}{ω^{' 2} - ω^{2}} d ω^{'} . \end{matrix}

(10)

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 (\cdot)

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

\dot{C}

S_{p} (L)

, overlayed on top of each other and sorted by the ascending index p, suspended above

\dot{C}

, and glued to each other along the slit extending from 0 to

- \infty

, which each sheet owns [23], as shown in Figure 1. If

z^{'}

is a point of

\dot{C}

and

\overset{˘}{z}

is the point belonging to

S_{q} (L)

, for the q-th sheet of

S (L)

, which lies above the point

z^{'}

\dot{C}

, the value taken by

L

\overset{˘}{z}

L_{\overset{˘}{z^{'}}}

, is given by [23]:

L_{\overset{˘}{z}} = l n | z^{'} | + i [a r g_{π} (z^{'}) + 2 π q] .

(11)

Based on the above, if

\overset{˘}{γ}

is the copy on

S (L)

of the path

γ (t)

lying in

\dot{C}

, we have from (11) that between the values taken by

L (\cdot)

z^{'} \in γ (t)

and the values assumed by

L

\overset{˘}{z} \in \overset{˘}{γ}

, the following relation must hold [23]:

L (z^{'}) = l n | z^{'} | + i [a r g_{π} (z^{'}) + 2 π p],

(12)

where, for each

\overset{˘}{z} \in \overset{˘}{γ}

, the value assigned to the integer parameter p is the index q of the Riemann sheet

S_{q} (L)

on which

\overset{˘}{z}

lies. Relation (12) can be used to compute

L (\cdot)

γ (t) ≜ S_{21} / (1 - S_{11} R)

. Substituting (8) into (12), we have the following:

Im [L (\frac{S_{21}}{1 - S_{11} R})] = [a r g_{π} (\frac{S_{21}}{1 - S_{11} R}) + 2 π q] .

(13)

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:

n_{e f f} k_{0} d = [a r g_{π} (\frac{S_{21}}{1 - S_{11} R}) + 2 π q]

(14)

Inserting (10) into (14), we obtain for the q values as follows:

q = \frac{1}{2 π} [k_{0} d (1 - \frac{2}{π} P \int_{0}^{+ \infty} \frac{ω^{'} κ_{e f f} (ω^{'})}{ω^{' 2} - ω^{2}} d ω^{'}) - a r g_{π} (\frac{S_{21}}{1 - S_{11} R})] .

(15)

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]:

\begin{matrix} \hat{q} & = \underset{q \in Z}{argmin} ϵ_{e r r} \end{matrix}

(16)

\begin{matrix} ϵ_{e r r} & = |a r g_{π} (\frac{S_{21}}{1 - S_{11} R}) + 2 π q - (1 - \frac{2}{π} P \int_{0}^{+ \infty} \frac{ω^{'} κ_{e f f} (ω^{'})}{ω^{' 2} - ω^{2}} d ω^{'}) k_{0} d| \end{matrix}

(17)

in which the q integers are those minimizing the error,

ϵ_{e r r}

, 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

a r g_{π} (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.

Algorithm 3 Numerical integration-like phase unwrapping

a : = 0

;

a r g_{α_{0}} (\cdot) : = a r g_{π} (\cdot)

;

γ_{t_{0}} : = \frac{S_{21} (a)}{1 - S_{11} (a) R (a)}

;

{[L (γ_{t})]}_{k = 0} = l n | γ (a) | + i a r g_{α_{0}} (γ (a))

;

5: for k: = 1 to n do;

t_{k} : = ω_{k}

;

γ_{t_{k}} : = \frac{S_{21} (t_{k})}{1 - S_{11} (t_{k}) R (t_{k})}

;

a r g_{α} : = Func_Kramers_Kronig_Integral

p_{k} : = Func_Minimize (a r g_{α}, a r g_{α_{0}} (γ_{t_{k}}))

;

10:

a r g_{α_{k}} (\cdot) : = a r g_{α_{0}} (\cdot) + 2 π p_{k}

;

11:

[L (γ_{t}))]_{k} : = l n | γ_{t_{k}} | + i a r g_{α_{k}} (γ_{t_{k}}))

;

12: end for

3. Numerical Results

To compare the performance of the phase-unwrapping approaches described in Section 2.2 and Section 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₃) 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₃ 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

a r g_{π} (\cdot)

, 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

a r g_{π} (\cdot)

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,

a r g_{π} (\cdot)

, 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).

4. 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₃ 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.

Author Contributions

Conceptualization, G.A. and M.V.; methodology, G.A.; software, G.A.; validation, S.C., M.V. and P.D.B.; formal analysis, G.A.; investigation, M.V.; resources, S.C.; data, P.D.B.; writing—original draft preparation, G.A., M.V. and S.C.; writing—review and editing, S.C.; visualization, S.C.; supervision, G.A. and M.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

$ϵ_{e f f}$	Effecive electric permittivity;
$μ_{e f f}$	Effective magnetic permeability;
$N_{e f f}$	Complex refractive index;
$n_{e f f}$	Effective refractive index;
$k_{0}$	Free space wavenumber;
d	MM thickness;
R	Reflection coefficient;
z	Effective impedance;
$S_{11}, S_{21}$	Scattering parameters;
$L O G (\cdot)$	Complex logarithm;
$l o g (\cdot)$	Principal logarithm function;
$l n (\cdot)$	Natural logarithm function;
$\| \cdot \|$	Absolute value function;
$a r g_{π} (\cdot)$	Principal argument function;
p	Branch index;
$e^{(\cdot)}$	Complex exponential function;
$L (\cdot)$	Right inverse function;
$L$	Global analytic logarithm;
$P$	Cauchy principal value;
$C$	Complex plane;
$\dot{C}$	Complex punctured plane;
$S (L)$	Riemann surface of $L$ .

Appendix A. Analytic Continuation, Phase Unwrapping, and Homomorphic Systems

In the homomorphic signal processing scope, Oppenheim and Schafer introduced phase unwrapping to make the homomorphic filtering operation possible [35]. Homomorphic filtering is founded on the concept of generalized linear superposition, which is defined as follows:

ϕ_{S} (c ⊙ ξ_{i} (τ) \oplus ξ_{j} (τ)) = c ⊡ ϕ_{S} (ξ_{i} (τ)) ⊞ ϕ_{S} (ξ_{j} (τ))

(A1)

where:

ϕ_{S} (\cdot) = (ϕ_{H^{″}} \circ ϕ_{H^{'}}) (\cdot)

(A2)

are proper input–output functions. The set

[⊙, \oplus]

is made up of an input generalized sum, which takes places between a pair of valid input elements

ξ_{i} (τ), ξ_{j} (τ)

, and by an input generalized multiplication, which is made between a scalar c and an input element

ξ_{i} (τ)

. The second set,

[⊡, ⊞]

, in a symmetrical way, is composed of an output-generalized sum and an output-generalized multiplication. They operate on scalars c and output function elements

ϕ_{S} (ξ_{i} (τ)), ϕ_{S} (ξ_{j} (τ))

. In the framework of the homomorphic filtering, the (A1) specializes as follows:

e^{(l n ({(ξ_{i} (τ))}^{c} \cdot l n (ξ_{j} (τ))))} = {(ξ_{i} (τ))}^{c} \cdot (ξ_{j} (τ))

(A3)

where

[\oplus, ⊙] = [⊡, ⊞] = [\cdot, {(\cdot)}^{c}]

and

ϕ_{H^{″}} = e^{(\cdot)}; ϕ_{H^{'}} = l n (\cdot)

. In this case, the uniqueness and the correctness of the input–output operations are guaranteed by the fact that

l n (\cdot)

is the inverse of

e^{(\cdot)}

. To guarantee the same in the case where

c \in C

and

ξ_{i} (τ)

ξ_{j} (τ)

are complex-valued functions, Oppenheim and Schafer developed a procedure to restore the continuity of the argument of the principal logarithm

l o g (\cdot) = l n | \cdot | + i [a r g_{π} (\cdot)]

, avoiding the jumps caused by the crosses with the

R^{-}

semi-axis, through the computation of the integer parameters p needed to compensate these discontinuities [26,35]. Considering the results provided by Theorem 1, we can understand that their procedure is exactly the analytic continuation defined by this theorem, which guarantees the uniqueness of the output operations of a multiplicative homomorphic system by ensuring that the function

ϕ_{H^{'}} (\cdot)

avoids the jumps of

[a r g_{π} (\cdot)]

caused by crosses with the

l o g (\cdot)

branch-cut, which is an admissible right-inverse function

l (\cdot)

for the complex exponential function

e^{(\cdot)}

References

Neelakanta, P.S. Handbook of Electromagnetic Materials: Monolithic and Composite Versions and Their Applications; CRC Press: Boca Raton, FL, USA, 1995. [Google Scholar]
Tong, X.C. Advanced Materials and Design for Electromagnetic Interference Shielding; CRC Press: Boca Raton, FL, USA, 2016. [Google Scholar]
Chen, L.F.; Ong, C.K.; Neo, C.P.; Varadan, V.V.; Varadan, V.K. Microwave Electronics: Measurement and Materials Characterization; John Wiley & Sons: Hoboken, NJ, USA, 2004. [Google Scholar]
Xiaojian, F.; Cui, T.-J. Recent progress on metamaterials: From effective medium model to real-time information processing system. Prog. Quantum Electron. 2019, 67, 100223. [Google Scholar]
Smith, D.R.; Padilla, W.J.; Vier, D.C.; Nemat-Nasser, S.C.; Schultz, S. Composite medium with simultaneously negative permeability and permittivity. Phys. Rev. Lett. 2000, 84, 4184. [Google Scholar] [CrossRef]
Smith, D.R.; Schultz, S.; Markoš, P.; Soukoulis, C.M. Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients. Phys. Rev. B 2002, 65, 195104. [Google Scholar] [CrossRef]
Vineetha, K.V.; Madhav, B.T.P.; Kumar, M.S.; Das, S.; Islam, T.; Alathbah, M. Development of Compact Bandpass Filter Using Symmetrical Metamaterial Structures for GPS, ISM, Wi-MAX, and WLAN Applications. Symmetry 2023, 15, 2058. [Google Scholar] [CrossRef]
Hossain, M.B.; Faruque, M.R.I.; Islam, M.T.; Khandaker, M.U.; Tamam, N.; Sulieman, A. Modified coptic cross shaped split-ring resonator based negative permittivity metamaterial for quad band satellite applications with high effective medium ratio. Materials 2022, 15, 3389. [Google Scholar] [CrossRef] [PubMed]
Feng, S.N.A.; Wang, Y.; Fei, S.; Yan, Z.; Yu, L.; Chen, J.; Liu, F. Dual ultrahigh-Q Fano Resonances of 3D gap metamaterials for slow light from ultraviolet to visible range. Opt. Commun. 2023, 549, 129811. [Google Scholar] [CrossRef]
Alam, M.J.; Latif, S.I. Double-Split Rectangular Dual-Ring DNG Metamaterial for 5G Millimeter Wave Applications. Electronics 2022, 12, 174. [Google Scholar] [CrossRef]
Pallavi, M.; Kumar, P.; Ali, T.; Shenoy, S.B. Modeling of a negative refractive index metamaterial unit-cell and array for aircraft surveillance applications. IEEE Access 2022, 10, 99790–99812. [Google Scholar] [CrossRef]
Angiulli, G.; Versaci, M.; Calcagno, S.; Di Barba, P. Quick retrieval of effective electromagnetic metamaterial parameters by using a Multi-fidelity Surrogate Modelling approach. Eur. Phys. J. Appl. Phys. 2020, 90, 20901. [Google Scholar] [CrossRef]
Gao, J.; Zhang, Y.; Wu, Q. All dielectric terahertz left-handed metamaterial based on mie resonance coupling effects. IEEE Access 2019, 7, 94882–94887. [Google Scholar]
Zhang, S.; Arya, R.K.; Whittow, W.G.; Cadman, D.; Mittra, R.; Vardaxoglou, J.C. Ultra-wideband flat metamaterial GRIN lenses assisted with additive manufacturing technique. IEEE Trans. Antennas Propag. 2020, 69, 3788–3799. [Google Scholar] [CrossRef]
Simovski, C.R. Material parameters of metamaterials (a review). Opt. Spectrosc. 2009, 107, 726–753. [Google Scholar] [CrossRef]
Simovski, C.R. On electromagnetic characterization and homogenization of nanostructured metamaterials. J. Opt. 2010, 13, 013001. [Google Scholar] [CrossRef]
Simovski, C.R. Bloch material parameters of magneto-dielectric metamaterials and the concept of Bloch lattices. Metamaterials 2007, 1, 62–80. [Google Scholar] [CrossRef]
Arslanagić, S.; Hansen, T.V.; Mortensen, N.A.; Gregersen, A.H.; Sigmund, O.; Ziolkowski, R.W.; Breinbjerg, O. A review of the scattering-parameter extraction method with clarification of ambiguity issues in relation to metamaterial homogenization. IEEE Antennas Propag. Mag. 2013, 55, 91–106. [Google Scholar] [CrossRef]
Cao, Z.X.; Yuan, F.G.; Li, L.H. An automated phase correction algorithm for retrieving permittivity and permeability of electromagnetic metamaterials. AIP Adv. 2014, 4, 067115. [Google Scholar] [CrossRef]
Shi, Y.; Li, Z.Y.; Li, L.; Liang, C.H. An electromagnetic parameters extraction method for metamaterials based on phase unwrapping technique. Waves Random Complex Media 2016, 26, 417. [Google Scholar] [CrossRef]
Singh, H.; Singh Sohi, B.; Gupta, A. Thickness invariant parameter retrieval techniques for permittivity and permeability measurement. J. Microw. Power Electromagn. Energy 2018, 52, 215. [Google Scholar] [CrossRef]
Angiulli, G.; Versaci, M. Retrieving the effective parameters of an electromagnetic metamaterial using the Nicolson–Ross–Weir method: An analytic continuation problem along the path determined by scattering parameters. IEEE Access 2021, 9, 77511. [Google Scholar] [CrossRef]
Angiulli, G.; Versaci, M. Extraction of the electromagnetic parameters of a metamaterial using the Nicolson–Ross–Weir method: An analysis based on global analytic functions and Riemann surfaces. Appl. Sci. 2022, 12, 11121. [Google Scholar] [CrossRef]
Rothwell, E.J.; Cloud, M.J. Electromagnetics; CRC Press: Boca Raton, FL, USA, 2018. [Google Scholar]
Markushevich, A.I. Theory of Functions of a Complex Variable, 2nd ed.; Chelsea Publishing Company: New York, NY, USA, 1985. [Google Scholar]
Mowlaee, P.; Kulmer, J.; Stahl, J.; Mayer, F. Single Channel Phase-Aware Signal Processing in Speech Communication: Theory and Practice; John Wiley & Sons: Hoboken, NJ, USA, 2016. [Google Scholar]
Kaplan, S.T.; Ulrych, T.J. Phase unwrapping: A review of methods and a novel technique. In Proceedings of the CSPG CSEG Convention, Calgary, AB, Canada, 14–16 May 2007; pp. 534–537. [Google Scholar]
Varadan, V.V.; Ro, R. Unique retrieval of complex permittivity and permeability of dispersive materials from reflection and transmitted fields by enforcing causality. IEEE Trans. Microw. Theory Tech. 2007, 55, 2224. [Google Scholar] [CrossRef]
Szabó, Z.; Park, G.; Hedge, R.; Li, S.L. A unique extraction of metamaterial parameters based on Kramers–Kronig relationship. IEEE Trans. Microw. Theory Tech. 2010, 58, 2646–2653. [Google Scholar] [CrossRef]
Yoo, S.; Lee, S.; Choe, J.H.; Park, Q.H. Causal homogenization of metamaterials. Nanophotonics 2019, 8, 1063–1069. [Google Scholar] [CrossRef]
Peiponen, K.E.; Vartiainen, E.M.; Asakura, T. Dispersion, Complex Analysis and Optical Spectroscopy: Classical Theory; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1998; Volume 147. [Google Scholar]
Peiponen, K.E.; Lucarini, V.; Vartiainen, E.M.; Saarinen, J.J. Kramers–Kronig relations and sum rules of negative refractive index media. Eur. Phys. J. Condens. Matter Complex Syst. 2004, 41, 61–65. [Google Scholar] [CrossRef]
Wheeler, M.S.; Aitchison, J.S.; Mojahedi, M. Three-dimensional array of dielectric spheres with an isotropic negative permeability at infrared frequencies. Phys. Rev. B 2005, 72, 193103. [Google Scholar] [CrossRef]
Ohta, K.; Ishida, H. Comparison among several numerical integration methods for Kramers–Kronig transformation. Appl. Spectrosc. 1988, 42, 952–957. [Google Scholar] [CrossRef]
Oppenheim, A.; Schafer, R.W. Discrete Time Signal Processing; Prentice-Hall: Kent, OH, USA, 2009. [Google Scholar]

Figure 1. A portion of the Riemann surface

S (L)

suspended on the unitary disk

D \in \dot{C}

Figure 1. A portion of the Riemann surface

S (L)

suspended on the unitary disk

D \in \dot{C}

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

r_{0}

, sphere radius).

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

r_{0}

, sphere radius).

Figure 3. Arrays of LiTaO₃ 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 3. Arrays of LiTaO₃ 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. Arrays of LiTaO₃ 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.

Figure 4. Arrays of LiTaO₃ 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.

Figure 5. Arraysof LiTaO₃ spheres,

d = 130 μ 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 5. Arraysof LiTaO₃ spheres,

d = 130 μ 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.

Table 1. CPU time (s),

d = 30 μ m

Table 1. CPU time (s),

d = 30 μ m

$n_{s}$	Alg1	Alg2	Alg3	K-K(B)
512	$0.329$	$0.329$	$0.349$	$0.297$

Table 2. CPU time (s),

d = 70 μ m

Table 2. CPU time (s),

d = 70 μ m

$n_{s}$	Alg1	Alg2	Alg3	K-K(B)
1024	$0.362$	$0.343$	$0.501$	$0.4981$

Table 3. CPU time (s),

d = 130 μ m

Table 3. CPU time (s),

d = 130 μ m

$n_{s}$	Alg1	Alg2	Alg3	K-K(B)
8192	$0.446$	$0.2936$	$4.027$	$3.955$

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Angiulli, G.; Versaci, M.; Calcagno, S.; Di Barba, P. Analytic Continuation, Phase Unwrapping, and Retrieval of the Refractive Index of Metamaterials from S-Parameters. Sensors 2024, 24, 912. https://doi.org/10.3390/s24030912

AMA Style

Angiulli G, Versaci M, Calcagno S, Di Barba P. Analytic Continuation, Phase Unwrapping, and Retrieval of the Refractive Index of Metamaterials from S-Parameters. Sensors. 2024; 24(3):912. https://doi.org/10.3390/s24030912

Chicago/Turabian Style

Angiulli, Giovanni, Mario Versaci, Salvatore Calcagno, and Paolo Di Barba. 2024. "Analytic Continuation, Phase Unwrapping, and Retrieval of the Refractive Index of Metamaterials from S-Parameters" Sensors 24, no. 3: 912. https://doi.org/10.3390/s24030912

APA Style

Angiulli, G., Versaci, M., Calcagno, S., & Di Barba, P. (2024). Analytic Continuation, Phase Unwrapping, and Retrieval of the Refractive Index of Metamaterials from S-Parameters. Sensors, 24(3), 912. https://doi.org/10.3390/s24030912

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Theory

2.1. The NRW Approach for Recovering the Complex Refractive Index

2.2. Analyticity and Phase Unwrapping

2.3. Riemann Surfaces and Phase Unwrapping

3. Numerical Results

4. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Abbreviations

Appendix A. Analytic Continuation, Phase Unwrapping, and Homomorphic Systems

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI